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I. INTRODUCTION 


Viscous incompressible flows have played an important role both in basic fluid dy- 
namics research and in industrial applications. In two-dimensional (2-D) studies, incom- 
pressible flows have been a convenient choice for laboratory experiments such as low speed 
wind-tunnel or water-tank experiments. Since computational study of these flow problems 
has been performed for several decades, various viscous incompressible flow solution meth- 
ods have been developed. In many, perhaps a majority of, realistic engineering applications, 
incompressible flows are also encountered. To name a few, there are problems related to 
low speed aerodynamics and hydrodynamics such as the flow around high lift devices, 
the flow around submerged vehicles, flow through impeller passages, mixing of the flow in 
chemical reactors, coolant flow in nuclear reactors, and certain biofluid problems. There 
are algorithmic simplifications as well as geometric modeling involved in computing these 
flows. Furthermore, significant physical modeling is also required, such as turbulence mod- 
eling for high-Reynolds number flows. Therefore, numerical computation of these flows, 
especially in three dimensions, becomes partly an art and partly a science, hence we use 
the term ‘numerical simulation’. 

To date, computational flow simulation has not been a critical element in resolv- 
ing many engineering problems. For instance, impellers, automobiles, submarines, and 
chemical reactors have been designed reasonably well by empirical means. As technologies 
advance, modern flow devices tend toward a more compact and highly efficient design. 
Therefore, the approach relying on empiricism and simplified analysis becomes inadequate 
for resolving problems associated with those devices thar require advanced analysis. For 
example, in analyzing and redesigning the current Space Shuttle main engine (SSME) 
power head, computational simulations became an economical and time-saving supple- 
ment to experimental data. There are vast numbers of other real-world problems which 
demand accurate viscous, incompressible flow solutions, such as rocket-engine fuel flow, 
flow through an impeller, and blood flow through a ventricular assist device. Therefore, 
it is of considerable interest to have a computational fluid dynamics (CFD) capability for 
simulating these applications as an alternative to analytical or empirical approaches. 

The present lecture note attempts to summarize recent progress in developing vis- 
cous incompressible flow solvers and their application to real world problems, both fun- 
damental and applied in nature. Because of the wide range of geometries of interest, 
flexibility and convenience becomes an important part of the flow solvers. Among various 
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approaches, the primitive variable formulation is considered to be the most promising for 
three-dimensional (3-D) applications and is therefore emphasized. In the present note, 
some details of derivation of the equations and algorithms are given as well as the physical 
interpretation of the solution procedures, mainly within a finite-difference framework. 

In Chapter II, a brief review on existing methods using primitive variables is given. 
In Chapter III, the pseudocompressibility method is explained in some detail followed by 
a description of a fractional step method in Chapter IV. Then, in Chapter V, various 
validation computations axe presented. Finally, in Chapter VI, examples of real world 
applications are presented. 
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II. SOLUTION METHODS 


In this chapter, a brief review is given on promising methods for general applications 
to 3-D real world problems. Naturally for 3-D applications, flexibility, robustness and 
computational efficiency as well as accuracy become very important. Other methods used 
in the past for more fundamental fluid dynamics problems can be found in the literature 
(see, for example, Roach, 1972; Orszag and Israeli, 1974). 


2.1 Formulation 


Unsteady, 3-D, incompressible flow with constant density is governed by the follow- 
ing Navier-Stokes equations: 


dui 


+ 


du{ 

dxi 

dUiUj 


= 0 


dp dri' 


dx i dxj 


( 2 . 1 ) 

( 2 . 2 ) 


dt ' dxj 

where t is the time, Xi the Cartesian coordinates, U{ the corresponding velocity components, 
p the pressure, and T,j the viscous-stress tensor. Here, all variables axe nondimensionalized 
by a reference velocity and length scale. The viscous stress tensor can be written as 


T\j — 2uS{j Rij 
•' “ 2' dx, + dxi’ 


(2.3) 

(2.4) 


Here, v is the kinematic viscosity, Sij is the strain-rate tensor, and Rij is the Reynolds 
stresses. Various levels of closure models for Rij are possible. In the present study, 
turbulence is simulated by an eddy viscosity model using a constitutive equation of the 
following form: 


Rij = ^ RkkSij - 2 u t Sij (2.5) 

where vt is the turbulent eddy viscosity. By including the normal stress, Rkk, in the 
pressure, v in equation (2.3) can be replaced by ( v + ut) as follows. 


Tij = 2(u + ut)Sij = 2u T Sij (2.6) 

In the remainder of this note the total viscosity, t*r, will be represented simply by v. 
Therefore, in describing solution methods and algorithms, the incompressible Navier-Stokes 
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equations are modified to allow variable viscosity in the present formulation. Hence, the 
flow solvers described here include turbulent flow cases as well. 


The major difference between the incompressible and the compressible Navier-Stokes 
formulation is the lack of a time derivative term in the continuity equation in the incom- 
pressible formulation. Therefore, satisfying the mass conservation equation is the primary 
issue in solving the above set of equations. Physically, incompressible flow is characterized 
by elliptic behavior of the pressure waves, the speed of which in a truly incompressible 
flow is infinite. The pressure field is desired as a part of the solution. However, the pres- 
sure cannot be obtained directly from the governing equations. One can utilize a stream 
function- vorticity formulation to eliminate the pressure. This method has been successfully 
used in many 2-D flow problems. Since an extension of this method to three dimensions is 
difficult, a vorticity- velocity method was developed to remove pressure from the solution 
procedure. This method eliminates the use of pressure, however, at the expense of intro- 
ducing vorticity boundary conditions which are the key to the development of the flow 
field and not straightforward to specify. In realistic 3-D problems, these derived quantities 
are difficult to define or impractical to use. The primitive variable formulation, namely, 
using pressure and velocities as dependent variables, then becomes more convenient and 
flexible in 3-D applications. However, in this formulation, mass conservation and its re- 
lation to pressure has to be properly handled while achieving computational efficiency. 
Various techniques have been developed in the past, none of which have been proven to be 
universally better than the other. 


2.2 Method Based on Pressure Iteration 
2.2.1 MAC Method 

Perhaps the first primitive variable method is the marker-and-cell (MAC) method 
developed by Haxlow and Welch (1965). In this method, the pressure is used as a mapping 
parameter to satisfy the continuity equation. By taking the divergence of the momentum 
equation, the Poisson equation for pressure is obtained: 

(2- 

where 

L _ chiiUj , drij 
* dxj + dxj 


V 2 p = 


dhi d aui 

dxi dt dxi 


= 9 
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The usual computational procedure involves choosing the pressure field at the cur- 
rent time step such that continuity is satisfied at the next time step. The original MAC 
method is based on a staggered arrangement on a 2-D Cartesian grid (figure 2.1). The 
staggered grid conserves mass, momentum, and kinetic energy in a natural way and avoids 
odd-even point decoupling of the pressure encountered in a regular grid (Gresho and Sani, 
1987). However, these properties become unclear when applied to generalized coordinates. 
Even though the original method used an explicit Euler solver, various time advancing 
schemes can be implemented here. Ever since its introduction, numerous variations of the 
MAC method have been devised and successful computations have been made. Many more 
examples can be found in the literature, for example, Roach (1972), Ferziger (1987), and 
Orszag and Israelii (1974). 

The major drawback of this method is the amount of large computing time required 
for solving the Poisson equation for pressure. When the physical problem requires a very 
small time step, the penalty paid for an iterative solution procedure for the pressure may 
be tolerable. But the method as a whole is slow and the pressure boundary condition is 
difficult to specify. 

One important aspect of the numerical solution of the Poisson equation for pressure 
is tied to the spacial differencing of the second derivatives. To illustrate this, equation 
(2.7) is solved in three different ways as follows: 

Method 1 : 

First, an exact form of the Laplacian operator is used in solving equation (2.7). The 
Fourier transform of equation (2.7) is 



9 ' 


(2.8 a) 


where 


p = Fourier transform of p 
k 2 = k x 2 + k v 2 + k 2 

k x ,k y , k z = wave numbers in the x—, y— , and z — directions 
g 1 = finite difference approximation to g 
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The wave number, fcj, in discrete Fourier expansion is defined as 


ki = _ _ , 7 i = wave number in the Xi — direction 

N A 

n = -N/ 2,J,l,..(JV/2-l) 

A = mesh space 

N = number of mesh points in the Xj — direction 
By inverse transformation, p can be obtained. 

Method 2: 

A second approach is to use the difference form of the second derivatives in equation (2.7), 

<£ + £ + £>'-*’ <»•»> 

The Fourier transform of the above equation is 

-kikip = g 

where ki is the Fourier transform of the difference form of the second derivative. 


Method 3: 

The finite difference form of the governing equations is 


8ui 

It 




( 2 . 1 ') 


= Dtij = 0 (2.2') 

0X{ 

where h\ is the finite difference form of h{. By applying the divergence operator, D, to the 
above equations, the following difference form of the Poisson equation is obtained. 


DGp = 


8Du{ 

It 


+ Dh\ = g'i 


Then, taking the Fourier transform of the above equation, 


—k'ik'ip = g' 


(2.8c) 


Since equation (2.1') is solved numerically, this equation is compared in the above 
three methods. Again, the Fourier transform of equation (2.1') is 
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t 


(2.9) 


Sin 

6i 



To satisfy the continuity equation in grid space, the following 



has to be satisfied for a flow that has Du » = 0 at the beginning. In Fourier space, this 
is equivalent to k[ui = 0 at the first and next time steps. Substituting p from the above 
three methods into equation (2.9), the following results are obtained. 


For Method 1: 

II 

*** 

<3 



For Method 2: 

km- 

(K- k -§-m*o 

(2.10) 

For Method 3: 

6 . 




Tt m)= 

(K - - = 0 



The error introduced by method 1 and method 2 can be seen by observing the magnitude 
of A? 2 , lb 2 and k' 2 . For the purpose of comparison, modified wave numbers, k and k 1 are 
given below when five-point fourth-order central differencing is applied. 


(fc*) 2 = [15 — 16coa(Afci) + cos(2A&t)] 

(fc/) 2 = ^ [65 — 16cos(Afcj) — 64cos(2Afc») + 16cos(2A&i) — cos(4Afcj)] 

I 


(2.11) 


The magnitude of these quantities are compared in figure 2.2. Method 3 satisfies the 
continuity equation at the next time step in grid space, and hence should be used for the 
pressure field solution. Therefore, in the solution method using a Poisson equation for 
pressure, the divergence gradient (DG) operator plays an important role in satisfying the 
mass conservation in grid space. As will be explained later, this strict requirement can be 
relaxed in other approaches. 


2.2.2 SIMPLE Method 

The major drawback of the MAC method is that a Poisson equation must be solved 
for pressure. A direct solver is practical only for simple 2-D cases. However, for 3-D 
problems, an iterative procedure is the best available choice. The strict requirement of 
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obtaining the correct pressure for a divergence free velocity field in each step significantly 
slows down the overall computational efficiency. Since, for a steady-state solution, the 
correct pressure field is desired only when the solution is converged, the iteration procedure 
for the pressure can be simplified such that it requires only a few iteration at each time step. 
The best known method using this approach is the Semi-Implicit Method for Pressure- 
Linked Equations (SIMPLE) developed by Caretto et al. (1972), see also Patankar and 
Spalding (1972) or Patankar (1980). 

The method begins with a guessed pressure p*, which is usually assumed to be 
p w at the beginning of the cycle. Then the momentum equation is solved to obtain an 
intermediate velocity tt,*. 


Ui* — Ui n = At[fcn(u n ,v n ,u*,v*) — ^-] 
The corrected pressure is obtained by setting 

P = P* + P' 

The velocity correction is introduced in a similar manner. 

m = Ui * + Ui 


( 2 . 12 ) 


(2.13) 


(2.14) 


Now the relation between the the pressure correction, p', and the velocity correction, u^, 
is obtained from a simplified momentum equation. First, the equation for p' and u\ is 
obtained from the linearized momentum equations. Then, by dropping all terms involving 
neighboring velocities, the following form of the pressure correction equation is obtained. 


u 


i 

i 



(2.15) 


where u; is a function of the particular differencing scheme chosen. Substituting equations 
(2.14) and (2.15) into the continuity equation, a pressure correction equation is obtained. 
This is equivalent to taking the divergence of equation (2.15). This procedure in essence 
results in a simplified Poisson equation for pressure, which can be solved iteratively line- 
by-line. 


The unique feature of this method is the simple way of estimating the velocity 
correction u\. This feature simplifies the computation but introduces empiricism into 
the method. Despite its empiricism, the method has been used successfully for many 
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computations. It is not the intention of this note to fully evaluate this method. Further 
details of SIMPLE, SIMPLER and other variations can be found in the literature, for 
example, see Patankar (1980). 


2.2.3 Fractional-Step Method 


The fractional-step method can be used for time-dependent computations of the 
incompressible Navier- Stokes equations (see Chorin, 1968; Yanenko, 1971; Marchuk, 1975). 
Here, the time evolution is approximated by several steps. Operator splitting can be 
accomplished in several ways by treating the momentum equation as a combination of 
convection, pressure, and viscous terms. The common application of this method is done 
in two steps. The first step is to solve for an auxilary velocity field using the momentum 
equation in which the pressure-gradient term can be computed from the pressure in the 
previous time step (Dwyer, 1986) or can be excluded entirely (Kim and Moin, 1985). In 
the second step, the pressure is computed which can map the auxiliary velocity onto a 
divergence-free velocity field. 


This procedure is illustrated by the following example of Mansour (private commu- 
nication, NASA Ames, 1986). 

Step 1: Calculate auxiliary or intermediate velocity, di, by 


^ /o r T n ff n— I ^ ^ yt 2/- i 7i 


(2.16) 


where 

H ‘ = 

The advection terms are advanced by a second-order Adams-Bashforth method while the 
viscous terms are modified to use an implicit approximate factorization scheme. Here 
the pressure gradient term is added to the procedure by Kim and Moin (1985), so as to 
minimize the pressure correction in each time step. 


Step 2: Solve for the pressure correction. 

The second step of the momentum equation can be written as 

t .n+l _ a . i r 

— = (<£ n+1 _ <**) 

At 2 Sxi KV * J 


Ui‘ 


where 


p n = 6 n - — V 2 <£ n 

P V 2 Re * 


(2.17) 
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This equation combined with continuity equation results in the following Poisson equation 
for the pressure correction. 


V 2 (<T +1 - 4> n ) 


2 S A 


At 6xi 


Ui 


(2.18) 


Once the pressure correction is computed, new pressure and velocities are calculated 
as follows: 


„n+l 


At 


= p " + (* n+1 - r ) - - <n 




(2.19) 


2 Sxi 


Others (for example, Dwyer, 1986) used essentially the same procedure illustrated 
above. One particular aspect of the fractional step method requiring special care is the 
intermediate boundary conditions. Orszag et al. (1986) discussed this extensively. As 
will be explained later in this note, Rosenfeld et al. (1988) devised a generalized scheme 
where physical boundary conditions can be used at intermediate steps. As with other 
pressure based methods, the efficiency of the fractional step method depends on the Poisson 
solver. A multigrid acceleration, which is physically consistent with the elliptic field, is 
one possible avenue to enhance the computational efficiency. The procedure described 
above was applied in Cartesian coordinates to calculate time dependent flows using either 
a staggered grid (Kim and Moin, 1985) or a regular grid (Dwyer, 1986). Extension of 
the staggered grid to generalized coordinates is not straightforward and involves features 
related to geometric conservation. More details of this extension will be described later. 


2.3 Pseudocompressibility Method 


Recent advances in the state of the art in CFD have been made in conjunction 
with compressible flow computations. Therefore, it is of significant interest to be able to 
use some of these compressible flow algorithms. To do this, the artificial compressibility 
method of Chorin (1967) can be used. In this formulation, the continuity equation is 
modified by adding a time-derivative of the pressure term resulting in 


1 dp du{ 
(3 dt + dxi 


( 2 . 20 ) 


where (3 is an artificial compressibility or a pseudocompressibility parameter. Together 
with the unsteady momentum equations, this forms a hyperbolic-parabolic type of time- 
dependent system of equations. Thus, fast implicit schemes developed for compressible 
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flows, such as the approximate-factorization scheme by Beam and Warming (1978), can be 
implemented. It is to be noted that the t no longer represents a true physical time in this 
formulation. Various applications which evolved from this concept have been reported for 
obtaining steady-state solutions (e.g., Steger and Kutler, 1977; Kwak et al. 1986, Chang 
et al., 1985, 1988a, b; Choi and Merkle, 1985). To obtain time- dependent solutions using 
this method, an iterative procedure can be applied in each physical time step such that 
the continuity equation is satisfied. Merkle and Athavale (1987) and Rogers and Kwak 
(1988, 1989) reported successful computations using this pseudo time iteration approach. 

In an incompressible flow, a disturbance in the pressure causes waves which travel 
with infinite speed. When pseudocompressibility is introduced, waves of finite speed result 
for which the magnitude of the speed depends upon the pseudocompressibility constant /?. 
In a true incompressible flow, the pressure field is affected instantaneously by a disturbance 
in the flow, but with pseudocompressibility, there will be a time lag between the flow 
disturbance and its effect on the pressure field. In viscous flows, the behavior of the 
boundary layer is very sensitive to the streamwise pressure gradient, especially when the 
boundary layer is separated. If separation is present, a pressure wave traveling with finite 
speed will cause a change in the local pressure gradient which will affect the location of 
the flow separation. This change in separated flow will feed back to the pressure field, 
possibly preventing convergence to a steady state. An extensive mathematical account of 
the pseudocompressibility approach is given by Temam (1979). 

This method has been extensively utilized in developing several incompressible 
Navier-Stokes codes at NASA Ames Research Center. A more detailed discussion of this 
method is given separately in the next chapter. 

2.4 Method Using Vorticity Transport Equations 

Among various other nonprimitive variable methods, the vorticity- velocity method 
is promising for 3-D applications. This method is briefly introduced here. 

A vorticity- velocity method was proposed by Fasel (1976) to study the boundary 
layer stability problem in two dimensions. A number of other authors used this approach in 
solving incompressible flow problems (i.e. Dennis et al. 1979; Gatski et al., 1982; Osswald 
et al., 1987; Hafez et al., 1988). However, a three-dimensional extension of this method 
has been limited to simple geometries to date. 


11 



Taking the curl of the above equation and using the continuity equation, the following 
Poisson equation for velocity is obtained. 

V 2 u = -Vxw (2.22) 

Equations (2.21) and (2.22) can be solved for the velocity and vorticity. Here, the pressure 
term and the continuity equation are removed at the expense of introducing three vorticity 
equations. This requires vorticity boundary conditions on a solid surface. 

As in the case of the fractional step approach, the computational efficiency of this 
approach depends on the Poisson solver. Overall performance in general 3-D applications 
remains to be seen. 
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III. PSEUDOCOMPRESSIBILITY METHOD 


As briefly introduced in the previous chapter, several incompressible Navier-Stokes 
solvers have been developed based on the pseudocompressibility method. In this chapter, 
physical characteristics of this method is examined followed by a more detailed discussion 
on solution procedures which were developed utilizing this method. 


3.1 Formulation and Its Physical Characteristics 

The use of the pseudocompressibility method as shown by equation (2.20) results in 
a system of hyperbolic-parabolic equations of motion. Physically, this means that waves of 
finite speed are introduced into the incompressible flow field as a medium to distribute the 
pressure. For a truly incompressible flow, the wave speed is infinite, whereas the speed of 
propagation of these pseudo waves depends on the magnitude of the pseudocompressibility. 
Ideally, the value of the pseudocompressibility is to be chosen as high as the particular 
choice of algorithm allows so that the incompressibility is recovered quickly. This has to 
be done without lessening the accuracy and the stability property of the numerical method 
implemented. On the other hand, if the pseudocompressibility is chosen such that these 
waves travel too slowly, then the variation of the pressure field accompanying these waves 
is very slow. This will interfere with the proper development of the viscous boundary layer, 
especially when the flow separates. For internal flow, the viscous effect is important for 
the entire flow field and the interaction between the pseudo pressure- waves and the viscous 
flow field becomes very important. Pseudocompressibility relaxes the strict requirement 
of satisfying mass conservation in each step. However, to utilize this convenient feature, 
it is essential to understand the nature of the pseudocompressibility both physically and 
mathematically. Chang and Kwak (1984) reported details of the pseudocompressibility, of 
which some key features are summarized below. 


For simplicity of analysis, the 1-D form of the governing equations are considered 


here. 


dp du 

at + % _0 


du du 2 

Incite 


dp 




dx 


(3.1a) 

(3.16) 


where r w is the normalized shear stress at the wall. The normal stress term in the stream- 
wise direction which contributes to the diffusion of the waves is neglected here. Neglecting 
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the shear stress term for simplicity, the following equation is derived from equations 3.1a 
and 3.1b. 


dp 


+ 


dp 


du 1 

dt (u ± c) dt ' Ydx ' (u ± c) dx 


du 1 

+ 


(u ± c) = 0 


The pseudo speed of sound, c, is found to be 

c = 

Relative to this sound speed, the pseudo Mach number, M, can be expressed as 

u 


u 

M = — = . 

c \/ti 2 + /? 


<1 


(3.2) 


(3.3) 


(3.4) 


It is clear that M is always less than 1 for all /? > 0. Therefore, pseudocompressibility does 
not introduce shock waves into the system and the flow remains subsonic with respect to 
the pseudo-sound speed. 

To study the main features of the wave propagation, the momentum equation is 
locally linearized. By cross differencing equations (3.1a) and (3.1b), the following equations 
are obtained: 


d 2 p 

dt* 

d 2 u 

dt* 


0 d 2 p d 2 p dr w 

+ 2u dtdx 13 dx * 13 dx 

(3.5a) 

d 2 u „d 2 u dr w 

(3.56) 

+ 2u dtdx P dx* ~ dt 

d as 


9 , t A 9 1 (p\ _ / Pdr^/fc \ 

m + ( c) dx J \u)-\-dr w /dt) 

(3.6) 


d d 

m +{u + c) te 

If the shear stress term on the right hand side of equation (3.6) were absent, characteristic 
equations for the linear waves would take a simple form expressed by 


£+(-+< 0 ^] (£)=<> 
s +( u - c 4 ] (*-) =0 


(3.7a) 

(3.76) 


Waves having a positive sign propagate downstream with a speed \u + c|, and waves having 
a negative sign travel against the stream with a speed |u — c|. 

The presence of the shear stress term, however, complicates the wave systems since 
the shear stress depends on the velocity field. The coupling between the pseudo pressure 
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waves and the vorticity spreading depends on their respective time scales. An order-of- 
magnitude analysis can be done here to obtain a general guideline for pseudocompressibil- 
ity. Suppose the distance from a point in the flow field such as a point on a flat plate or 
channel to the downstream boundary is xl, then the time required for the upstream prop- 
agating wave to reach that point from the down stream boundary, t^, can be estimated 
by the following relation: 


XL-~ 



%ref 



U 


!L\ 

2 r L 

ref 


(3.8) 


On the other hand, the time scale, Tg, for the viscous effect to spread for a distance, 8, can 
be estimated as follows: 


£ „ 2 Vvi _ 2 v fare/ _ 2 jjs_ 

%ref %ref y u ref^ref ®rc/ V 


(3.9) 


In order for the viscous boundary layer to adjust to its new pressure environment properly 
and to avoid slow fluctuations of the separation region when it is present, it is required 
that 


Tg > Tt 


which leads to 




Sr ef 

8 



Xref 



(3.10) 


This gives an estimate of the lower bound of the pseudocompressibility parameter, /?. The 
physical phenomena described above will be illustrated by numerical experiments later in 
this report. 


3.2 Governing Equations in Generalized Coordinates 

To perform calculations on 3-D arbitrarily shaped geometries, the following gener- 
alized independent variables are introduced which transform the physical coordinates into 
general curvilinear coordinates 

£ = £(x,y,z,t) 

v = (3.H) 
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The Jacobian of the transformation is .defined as 


where 


J - 


d{x,y,z) 


£x Cy 

Vz Vy Vz 

C* Cv Cz 


u = 


d( _ dji 

dx ’ dy 


In actual coding, metric terms Eire calculated as follows : 


(U \ if Vv z < - y< z v \ 

I £v ] = j7 I x < z v - x v z < I 
\tzj \x v y ( -x < y v J 


fVx\ 1 ( VCH ~ ViH \ 
% = 77 x * z < ” *< z € ) 

\vzj - x m / 


f (*\ l f V( z v - VvH \ 

<v I = j7 \ x nH ~ x ( z v J 
\CzJ \xtyr, - z v y ( ] 


J' = 


d(s,y,z) 


x t x v x ( 

Vi y v y< 

Z( z v z ( 


Applying the transformation to equation (2.1) and (2.2) yields 



-% (e - ^ • t, (} ~ / ” ) - wfi - M = - f 


d 

d£ 



(3.12) 


(3.13) 


(3.14) 


(3.15) 


(3.16) 
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where 


„ 1 
u = — 

J 


. 1 
e== 7 


1-1 


. i 

* = 7 


V 

1U 


£®p + uU 
(yP + VU 
[tzP + wU 

rixp + uV 
r]yp + vV 
[r} z p + wV m 

C*P + uW 
CyP+vW 

ltzP + wW 


V = Tf t + r} x u + rjyV + 7} z w 
W = Ct + C*« + CyV + ( Z W 


( 3 . 17 ) 


The viscous terms axe quite lengthy and therefore given here separately. 


d 


d 


-Tij — 


dxj 1 dxj 

_ d_ 

dx 


uSi 


U x +U Z 

d 

Uy 

+ V X ' 

d 

u z + w x ' 

V z +Uy 


Vy 

+ Vy 

+ d~z 

V z + Wy 

w x + u Zm 

Wy 

+ »*. 

w z + w z 


where 


du du 

Ux = d^ ,Uv = d^ 


When v is constant, the contribution of the second term in each bracket sums up to be 
zero when the velocity field is divergence free. However, since v varies in space and time in 
general, theses terms have to be kept. Then the viscous terms in transformed coordinates 
are given by 


e. = ^vnv4 + v4 + v4) 




(V,-(V t J + V,J + V<£) 


u 

V 

w 


u 

V 

w 


/j. du dv dw x 

+ * m Wi *'Wi+ u Wi } 


, du dv dw v 

+ W'du +v, Wi +v ‘Wi’ 


f 

hi jj 

T hi 

d t, 

f: i 
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*--<vc.(v«| + v,J + vc£) 


u 

V 

w 


I. du dv dw^ 
+ (Cx^7+Cy^7 + Cz^:j 



The above full viscous terms can be simplified under special conditions. If an orthogonal 
coordinate system is used, then 


V&.Vfc = 0 ; for i ^ j 


For constant v (i.e., laminar flow), the contribution by the second term in parentheses in 
the above equation will approach zero as the steady state is reached, where equation (2.1) 
is numerically satisfied. Therefore, for flow with constant v in orthogonal coordinates, 
equation (3.18a) can be reduced as follows: 


Cv — (-j) (£as 2 + £y 2 + fz 2 ) 
h = (j)(Vx 2 +Vy 2 +Vz 2 ) 
i7, = (y)(<, 2 + <v 2 + <z 2 ) 


U( 

V( 

W£ 



u < 

v < 

w < 


(3.186) 


3.3 Steady-State Formulation 


The pseudocompressibility relation is introduced by adding a time derivative of 
pressure to the continuity equation, resulting in 


1 dp dui 

0 ^ + aTi 


(3.19) 


In the steady-state formulation the equations are to be marched in a time-like fashion until 
the divergence of velocity in equation (3.19) converges to zero. The time variable for this 
process no longer represents physical time, so in the momentum equations t is replaced 
with r, which can be thought of as a pseudo-time or iteration parameter. Combining 
equation (3.19) with the momentum equations gives the following system of equations: 

±D = -|(i - E.) - £(* - K) - |(G - 6.) = -R (3.20) 
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where R is the right-hand-side of the momentum equation and can be defined as the 
residual for the steady-state computations, and where 




P' 

U 

V 

Lt». 


E = 


F = 


G = 

E v = 
F v = 
G v = 


'0(U-t*)/J 1 = 1 
e J 


’Piy-ruyj l i 

/ J J 


'/3(W-Ct)/J ] 1 

9 J J 

' 0 ' 

' 0 ' 

M 

' o| 

9 v 


ix? + uu 
lyP + VU 
.£zP + wU. 
■0(V-rj t y 
7)xP + uV 
VyP + vV 

L t\ z p + wV . 

■m-Ct)- 

CxP + uW 
CyP + vW 

L CzP + vjW . 


For flow with constant v in orthogonal coordinates, 


where 


Ev = ( j) ft.* + £„ 2 + = Tl B 

Ev = ( j) {vx 2 + Vy 2 + Vz 2 )lm^ = 72 D 
^=(j)(C. 3 + Cv 2 + C 2 )^=73l? 


■o 

0 

0 

0- 

0 

1 

0 

0 

0 

0 

1 

0 

.0 

0 

0 

1. 


(3.21) 


(3.18c) 
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3.4 Steady-State Algorithm Using Approximate Factorization 


3.4.1 Difference Equations 

There are a number of different ways of defining variables in a grid system. Either standard 
cell-node oriented or a staggered arrangement can be chosen. In Cartesian coordinates, a 
staggered grid has several favorable properties. In generalized coordinates, these advan- 
tages become obscured because of the interpolation required. However, a fully conserva- 
tive differencing scheme can be devised which maintains the convenience of a staggered 
arrangement in a Poisson solver (Rosenfeld et al., 1988). Using any grid system, spatial 
differencing can be done either in finite-difference (Steger and Kutler, 1977) or finite- 
volume form (Yoon and Kwak, 1988). The finite- volume scheme usually behaves better 
near geometric singularities. In spacial differencing, both central differencing and upwind 
differencing have been tried. In this section, the work done using centred differencing will 
be discussed first. Most of the results presented later in this report were obtained using 
finite-difference schemes. This chapter will focus on algorithms based on a finite-difference 
approach. 


The pseudo-time derivative is replaced by a trapezoidal rule finite-differencing 
scheme resulting in 


D n+1 =D n + 


At 


\ot J \ or 


n+li 


+ 0(Ar 3 ) 


(3.22) 


where the superscript n refers to the n th pseudo-time iteration level. By substituting 
equation (3.20) into equation (3.22), one obtains 

D n+1 + ^ J[6t(E - E v ) n+1 + S V (F - F v ) n+1 + 6 ( {G - G v ) n+1 ] 

2 a (3.23) 

= D n - =^J[8((E - E v ) n + 6 V (F - F v ) n + 8 C (G - G v ) n ] 

The problem is to solve for D n+1 , and this is nonlinear in nature since E n+1 = E(D n+1 ) is 
a nonlinear function of D n+1 as are F n+1 and G n+1 . The following linearization procedure 
is used. A local Taylor expansion about u n yields 

E n+ 1 = E n + A n (D n+1 - D n ) + 0(Ar 2 ) 

pn+1 = pn + B n (D n+ 1 - D n ) + 0(Ar 2 ) (3.24) 

G n+ 1 = G n + C n (D n+1 - D n ) + 0(Ar 2 ) 
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where A, B, and C are the Jacobian matrices 


■ 0 

L x /3 

L 2 0 

Li0 ' 

Lx 

Q + Lxu 

L 2 u 

L 3 U 

L 2 

L\v 

Q + L 2 v 

L 3V 

L 3 

L\w 

L 2 w 

Q + L3W. 


dE * dF ~ dG , 

A =m’ B =W' c =m (3 - 26) 

The Jacobian matrices can all be represented by the following, 


1 x I ±J \ ^ T^IU- ^ 2 <* ^ 3 “ I /o nM 

Ai — j I r.-« n _l r.« I {o-to) 


A A A A 

where A, = A, £, or <7 for z = 1, 2, or 3, respectively, and 

Q = Lo + L\u + L 2 v + L3W 
L 0 = «i)«, £. = «i).. £2 = (fi),. £3 = (&), 

& = ((,<)■ oiC)for(i,B, or(?) 

Substituting equation (3.24) into equation (3.23) results in the governing equation in delta 
form 


{ / + 1 j [« e (i n - r,) + s v {B n - r 2 ) + s c (d n - r.)] | (z> n+1 - D n ) 

k(JB - Ev) n + S V (F - F v ) n + 8 {(G - (?t,) n l 


= -ArJ 


(3.27) 


where 


r ‘ D ” +1 = (j)v« • (vf | + v,£ + 


' 'dri 

h = At for trapezoidal differencing 
= 2 At for Euler implicit scheme 

At this point it should be noted that the notation of the form [£((A — T)]Z) refers to 

BA 

K 


j- ( (AD) - l ( (TD) and not $ fD - f D. 


3.4.2 Approximate Factorization 
ADI scheme: 

The solution of equation (3.27) would involve a formidable matrix inversion problem. 
With the use of an alternating direction implicit (ADI) type scheme, the problem could be 
reduced to the inversion of three matrices of small bandwidth, for which there exist some 
efficient solution algorithms. The particular ADI form used here is known as approximate 
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factorization (AF) (Beam and Warming, 1978). It is difficult to apply the AF scheme to 
equation (3.27) in its full matrix form. Noting that at the steady-state the left-hand side 
of equation (3.27) approaches zero, a simplified expression for the viscous term as shown 
in equation (3.18b) is used on the left-hand side. To maintain the accuracy of the solution, 
the entire viscous term is used on the right-hand side. Using these terms, the governing 
equation becomes 

I«L,£ < (0 w+1 - D n ) = RHS (3.28a) 

where 

L <=[ 1+ 

L, = (3.286) 

L<= [ I +l Ijn+1<<(<5n “ 73) 

and RHS is the right-hand side of equation (3.27). When second-order centred differencing 
is used, the solution to this problem becomes the inversion of three block tridiagonal 
matrices. The inversion problem is reduced to the three inversions 

(. L V )AD = RHS 

{L ( )AD = AD (3.29) 

(L()AD n+1 = AD 

These inversions are carried out for all interior points, and the boundary conditions can 
be implemented explicitly. It is possible, however, to implement the boundary conditions 
implicitly, which will be discussed in a later section. 

A guideline for estimating the lower bound of ft was given by equation (3.10) which 
was derived from physical reasoning. Following the guideline, it is desirable to choose (3 
as large as possible to make the pressure wave faster. There is, however, a bound of /? 
which comes from the particular algorithm chosen here, namely, the error introduced by 
the approximate factorization. In implementing the AF scheme leading to equation (3.28), 
the following second-order cross product terms axe introduced into the equation: 

This term must be kept smaller than the original terms in the equation. Including only 
the terms which contain /?, this restriction can be expressed as 


22 


or 


< l 


Recalling the expression for A” given by equation (3.26), the terms which have /? in them 
give the following 


At 






<1 


The term to the right of (3 in this inequality is essentially the change in in either 
the £, 77 , or £ direction. An estimate of the order of magnitude of this term for the grids 
generally used is given by 

which puts the restriction on (3 

0((3At) < 1 (3.30) 

For most problems, the restrictions for /? given by equations (3.10) and (3.30) are satisfied 
with a value for j3 in the range from 1 to 10. As will be shown later, this restriction on 
the upper value of /? can be relaxed if the factorization error involving j3 is removed, for 
example, by a line relaxation scheme. 

In applying the present approximate factorization scheme, it has been found that 
the stability of the scheme is dependent on the use of some higher-order smoothing terms. 
These are used to damp out higher frequency oscillations which arise in the solution because 
of the use of second-order central differencing, and will be discussed later. 


LU-SGS scheme; 


Recently, Yoon and Jameson (1987) developed an implicit lower-upper symmetric- 
Gauss-Seidel (LU-SGS) scheme for the compressible Euler and Navier-Stokes equations. A 
similar scheme is devised for the pseudocompressible formulation (Yoon and Kwak, 1989). 
This LU-SGS scheme is not only unconditionally stable but completely vectorizable in 
three dimensions. Spacial differencing can be either central or upwind depending on the 
numerical dissipation model which augments the finite volume method (Yoon and Kwak, 
1988). This scheme is briefly described below. 


Starting from an unfactored implicit scheme similar to equation (3.27) 


Hi 


S ( A + S V B + 8<C 


] j (D n+1 - D n 


) 


= - At 


6i(E - E v ) + 8 n (F - F v ) + S^G - G v ) 


(3.27') 
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(3.31a) 


The LU-SGS implicit factorization scheme can be derived as 

LtL^Lu = RHS 

where 

Li = I + ^{S~ ( A + + S~ V B + + 6~ ( C + — A~ — B~ — C~) 

A 

L d = I+ £(i + - A+ + B* - B- + <7+ - C-) (3.314) 

A 

L U =I+ \{S + tA~ + 6+ v B~ + S + ( C~ + A + + B+ + C + ) 

where S~( and are the backward- and forward-difference operators respectively. In the 
framework of the LU-SGS algorithm, a variety of schemes can be developed by different 
choices of numerical dissipation models and Jacobian matrices of the flux vectors. Further 
details of this method will be given in a later report (Yoon and Kwak, 1989). 

3.4.3 Numerical Dissipation or Smoothing 

Higher-order smoothing terms are required to make the present algorithm stable. 
These terms help to damp out the higher-order oscillations and odd- and even-point de- 
coupling in the solution which are caused by the use of central differencing. The smoothing 
term can be related to an upwind finite-difference approximation. The idea of splitting 
the upwinding scheme into the central differencing plus dissipation was successfully im- 
plemented into the first- and the third-order dissipation model by Jameson et al. (1981). 
Pulliam (1986) discussed an implicit dissipation model extensively. More recently Yoon 
and Kwak (1988) unified the dissipation models in the framework of a finite volume total 
variation diminishing (TVD) method. These models are being extended into the pseudo- 
compressible formulation (Yoon and Kwak, 1989). Here, only those specifics are discussed 
that are relevant to the constant coefficient model, which was originally used in conjunction 
with the pseudocompressible formulation. 

By including these smoothing terms, equations (3.28a) and (3.28b) become 

L(L v L((D n+1 — D n ) = RHSof (3.27) - e e [(V*A*) 2 + (V„A n ) 2 + (V c A c ) 2 ]D n (3.28) 
where 

— J+ — J n+ 1 8 t(A* - 71) + 

L v = J+^J n+1 ^(i?-7 2 ) + eiV,A n 
L<= /+^J n+1 *<(i?-73) + *iV c A< 
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Here, V and A represent forward and backward spacial-differencing operators, respectively. 
To preserve the tridiagonal nature of the system, only second-order smoothing can be used 
on the left-hand side of the equation, whereas fourth-order smoothing is used on the right- 
hand side. When the diagonal algorithm (described in a later section) is used, however, 
it is feasible to increase the bandwidth of the system to a pent adi agonal. This makes it 
possible to use fourth-order smoothing on the left-hand side of the equations also. The AF 
algorithm will be stable if e, and e e satisfy a certain relation (see Pulliam, 1986; Jameson 
and Yoon, 1986) as discussed below. 

To study the nature of the numerical smoothing, a 1-D form of the dissipation terms is 
represented as follows: 


[1 - £i V { A { ] (p" +1 - p n ) = - e .(V { A;)V (3.32) 

Suppose p is represented by the discrete Fourier expansion 

P = '£Hky* (3-33) 

n 

where 

p = Fourier transform of p 

k = —n = wavenumber 
N A£ 

n = —N/2, ...0, 1, ..(JV/2 — 1) 

N = number of mesh points 
Equation (3.32) can be written as 


[1 - £i fc'] - p") = -MV)’#” 


where 


k' = —2 + 2 cos(k) 

(k 1 ) 2 = 6 — 8cos(fc) + 2cos(2k) 
From this, the amplification factor can be defined as 


(3.34) 


[i - £i fc' - £ ,(t’) 2 ] 

” P n [1 - £ (i'] 

To damp out the numerical fluctuations as time advances, the absolute value of the am- 
plification factor a has to be less than one for all possible frequencies, i.e., 



<r| < 1 
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Noting that k' is always negative, this requirement leads to the following relation. 

e e < 2(1 - €ik') (3.36a) 

It can be shown that the above inequality is always satisfied if 

2e e < e; (3.366) 

The exact relation between these two coefficients can be determined only by a nonlinear 
stability analysis. In the present study, e; is taken to be three times larger than e e . From 
the expression given in equation (3.35), it is clear that if e* is too large, the rate of damping 
will be diminished. It may not be advantageous to take a very large value for ej over c e . 
The choice of e e depends on the Reynolds number and the grid spacings. However, as 
discussed later, large values of e t adversely affect the accuracy of the continuity equation. 
Therefore, the magnitude of c e is usually taken to be small. If grid sizes axe fine enough 
to resolve the changes in the flow field, it can be taken as small as 10 -3 . 

There are two major sources of inaccuracy associated with the numerical dissipation terms; 
namely, (1) the numerical dissipation terms effectively change the Reynolds number of 
the flow, and (2) the explicit smoothing terms added to the continuity equation do not 
conserve mass. In particular, the explicit smoothing on the pressure can affect whether 
or not the solution converges to an incompressible flow field. Chang and Kwak (1984) 
showed that the pseudo-pressure waves decay exponentially with time, and vanish as the 
solution converges. Thus the change in pressure with time approaches zero. When there is 
no explicit smoothing added to the equation, the divergence of the velocity field identically 
approaches zero. However, when explicit smoothing is included, as the change in pressure 
approaches zero, the divergence of the velocity approaches 

S: ■* ~jb A <> 2 + + < v < A () J ] t ( 3 - 37 > 

where is a difference form of the divergence operator and c el is the explicit smoothing 
parameter for the pressure. If e e is scaled by h, i.e., c e = A re e , equation (3.37) becomes 
independent of the time step. Depending on the magnitude of /3 and the local pressure 
gradient, this term can deteriorate the conservation of mass. However, it is to be noted 
that, in the numerical computations, the central differencing scheme is modified to include 
numerical dissipation terms resulting in what is essentially an upwinding scheme. There- 
fore, in the present approach, the mass conservation is satisfied in an upwinding sense. In 
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numerical experiments, we have seen that the divergence of the velocity field evaluated 
using second order differences does not decrease as fast as the residual terms which can be 
regarded as the divergence of the velocity in an upwinding sense. 

3.4.4 Diagonal Algorithm 

In a diagonal algorithm, a similarity transform is used to diagonalize the Jacobian 
matrices and uncouple the set of equations. The equations can then be solved by solv- 
ing scalar tridiagonal matrices instead of solving block tridiagonal matrices. A similarity 
transform which symmetrizes and diagonalizes the matrices of the compressible gas dy- 
namic equations has been used by Warming et al. (1975) and Turkel (1973). This method 
was used by Pulliam and Chaussee (1981) to produce a diagonal algorithm for the Euler 
equations. This method can be applied to the compressible Navier-Stokes equations to 
obtain a considerable savings in computing time (Flores, 1985). In this section, similar- 
ity transforms for the matrices used in the method of pseudocompressibility are presented. 
They are used in a diagonal algorithm, which results in a substantial reduction in computer 
time (Rogers et al., 1987). 

Similarity transformations exist which diagonalize the Jacobian matrices 

Ai = TikiTr 1 (3.38) 

where A* is a diagonal matrix whose elements are the eigenvalues of the Jacobian matrices 
and which is given by 

Q 0 0 0 

7 _ 0 Q 0 0 

{ 0 0 Q-Lo/2 + c 0 

. 0 0 0 Q - To/2 - c 

and where c is the pseudospeed of sound, which is given by 

c = yAq + Io/2) 2 +£(£?+ 1! +I|) (3.40) 

The Ti matrix is composed of the eigenvectors of the Jacobian matrix. For more details 
on the derivation of Ti, its inverse and eigenvectors, see Rogers et al. (1987). 

The implementation of the diagonal scheme involves replacing the Jacobi an matrices 
in the implicit operators with the product of the similarity transform matrices and the 
diagonal matrix as given in equation (3.38). The identity matrix in the implicit operators 


(3.39) 
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is replaced by the product of the similarity transform matrix and its inverse. A modification 
is made to the implicit viscous terms by replacing the I m matrix with an identity matrix so 
that the transformation matrices may also be factored out of these terms. This implicitly 
adds an additional viscous dissipation term to the pressure. The transformation matrices 
are now factored out of the implicit operators to give 

c, = r,[j+ ^ j<„( a, - -Wir- 1 (3.4i) 

where the implicit viscous terms are now given by 

7* = • VfilSt, (3.42) 

Since the transformation matrices Eire dependent on the metric quantities, factoring 
them outside the difference operators introduces an error. No modification has been made 
to the right-hand side of the equation and therefore these linearization errors will not affect 
the steady-state solution, only the convergence path toward the solution. 

The implementation of this algorithm over the block algorithm will result in a sub- 
stantial reduction in computational time per iteration because of the decrease in the num- 
ber of operations performed. Additionally, considerably less memory is required to store 
the elements on the left-hand side. This additional memory was used to further vectorize 
the existing code as follows. Since the solution of a tridiagonal block or scalar matrix is 
recursive, it is not vectorizable for loops which use the current sweep direction as the inner 
do-loop index. However, if a large number of these matrices are passed into the inversion 
routines at once, then vectorization can take place in the ‘nonsweep’ direction. 

3.4.5 Boundary Conditions 

Once the flow solver is developed, a boundary condition procedure has to be devised 
to be compatible with the solution algorithm. Boundary conditions play an important part 
in determining the overall accuracy, the stability property, and the convergence speed of 
the solution. There are different types of boundaries encountered in numerical simulations 
which include solid surface, inflow and outflow, and far-field boundaries. 

Solid Surface : 
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At a solid surface boundary, the usual no-slip condition is applied. In general the 
grid point adjacent to the surface will be sufficiently close so that constant pressure normal 
to the surface in the viscous boundary layer can be assumed. For a £ = constant surface, 
this can be expressed as 

(!L =o - (343) 

This boundary condition can be implemented either explicitly or implicitly. The implicit 
implementation, however, will enhance the stability of the code. This can be done during 
the £ sweep by including the following in the matrix to be inverted. 

I4fl w + £ADj,n,2 = f (3.44) 


where 



•-1 

0 

0 

0- 


'PL =2 — PL= r 


0 

0 

0 

0 


0 

c = 

0 

0 

0 

0 

/ = 

0 


. 0 

0 

0 

0. 


0 


Inflow, Outflow and Far-field Conditions : 


The inflow and outflow boundary conditions for an internal flow problem and the 
far-field boundary conditions for an external flow problem can be handled in much the 
same way. The incoming flow for both problems can be set to an appropriate constant 
as dictated by the problem. For example, at the inlet to a pipe, the pressure can be set 
to a constant and the velocity profile can be specified to be uniform. The downstream 
conditions, however, are the most difficult to provide. Simple upwind extrapolation is not 
well-posed. The best convergence rate is obtained if global mass is conserved. So to ensure 
the best results, the velocities and pressure are first updated using a second-order upwind 
extrapolation. For an exit at L = LMAX this is written 



Q 


n+1 

Imox-1 


( 4 * 2 ^ _ 0 n+1 


A Z2 

Azi 


-1 


(3.45) 


where 


Azi 


Zlmax 1 


A 22 


Zlmaz Zlmaz — 2 


Then, these extrapolated velocities are integrated over the exit boundary to obtain the 
outlet mass flux 


fh ou t 



(3.46) 
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Then the velocity components are weighted by the mass flux ratio to conserve global mass 


yn + 1 _ yh 

™ 0 ut 


(3.47) 


If nothing further is done to update the boundary pressure, this cam lead to discontinuities 
in the pressure because momentum is not being conserved. A method of weighting the 
pressure by a momentum correction was presented by Chamg et al. (1985), where the 
pressure condition is obtained by the mass weighted velocities 


,n+1 = ?h - z [ {wW)n+1 - {wW)i "1 + £ (vc ' v ° (If) +1 ' ( 


dwy 

d<) 


(3.48) 


where W is the contravariamt velocity. In obtaining this formula, it has been assumed 
that the streamlines near the exit plane are nearly straight. Any appreciable deviation 
will cause a discontinuity in the pressure and may lead to an instability. To avoid this, a 
momentum weighted pressure was used. This was obtained by integrating the momentum 
corrected pressure p n+1 and the extrapolated pressure p ft across the exit 



The final outlet pressure is then taken to be 



(3.49) 


(3.50) 


Under these downstream boundary conditions, global conservation of mass and momentum 
are ensured. Many practical applications have been solved using the current procedure. , 
However, the nonreflecting type boundary conditions according to Rudy and Skrikwerda 
(1980) may enhance the convergence speed. 


3.5 Time- Accurate Formulation 

Time dependent calculation of incompressible flows are especially time consuming 
due to the elliptic nature of the governing equations. Physically, this means that any local 
change in the flow has to be felt by the entire flow field. Numerically, this means that in each 
t.im* step, the pressure field has to go through one complete steady-state iteration cycle, 
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for example, by Poisson-solver-type pressure iteration or pseudocompressibility iteration 
method. In transient flow, the physical time step has to be small and consequently the 
change in the flow field may be small. In this situation, the number of iterations in each 
time step for getting a divergence-free flow field may not be as high as regular steady-state 
computations. However, the time-accurate computations are in general extremely time- 
consuming. Therefore, it is particularly desirable to develop computationally efficient 
methods either by implementing a fast algorithm and by utilizing computer characteristics 
such as vectorization and parallel processing. In this section, a time-accurate method using 
pseudocompressibility developed by Rogers and Kwak (1988, 1989) is described. 


In this formulation the time derivatives in the momentum equations are differenced 
using a second-order, three-point, backward-difference formula 


3tt n+1 _ 4£ n + 

2A t 


M-l+1 


(3.51) 


where the superscript n denotes the quantities at time t = nAf and f is the right-hand 
side given in equation (3.15). To solve equation (3.51) for a divergence free velocity at the 
n+1 time level, a pseudo-time level is introduced and is denoted by a superscript m. The 
equations are iteratively solved such that u n+1,m+1 approaches the new velocity u n+1 as 
the divergence of •u n+1,77l+1 approaches zero. To drive the divergence of this velocity to 
zero, the following artificial compressibility relation is introduced: 


pn+l,m+l _pn+l,m 

At 


— -fiV ■u n+1,m+1 


(3.52) 


where r denotes pseudo-time and f3 is an artificial compressibility parameter. Combining 
equation (3.52) with the momentum equations gives 


I tr (£) n + 

= _£n+l,m+l _ - 2 D n + 0.5 D n_1 ) 

At 


(3.53) 


where D is the same vector defined in equation (3.21), R is the same residual vector defined 
in equation (3.20), It T is a diagonal matrix, and I m is a modified identity matrix given by 


I tr = diag 


'J_ L5 U5 O 
At’ At’ At’ At 
I m = diag[0, 1, 1, 1] 
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Finally, the residual term at the m+1 pseudo- time level is linearized giving the following 
equation in delta form 


Itr (dR 

J + \dD 


n+l»ni' 




(3.54) 


-£ n+1 ’ m _ _^(i ,5D n+hm - 2 D n + 0.5 2> n-1 ) 
At ' 


As can be seen, this equation is very similar to the steady-state formulation given by 
equation (3.27), which can be rewritten for the Euler implicit case as 


1 

JAt 


1 + 



(D n+1 - D n ) = -R n 


(3.55) 


Both systems of equations will require the discretization of the same residual vector R. 
The derivatives of the viscous fluxes in this vector are approximated using second-order 
central differences. The convective flux terms can be discretized using central differences 
as was done in section 3.4. This will require numerical dissipation terms for stability. 
Since, in pseudo compressible formulation, the governing equations are changed into the 
hyperbolic-parabolic type, sdme of the upwind differencing schemes which have recently 
been developed for the compressible Euler and Navier-Stokes equations by a number of 
authors (i.e., Roe, 1981; Chakravarthy and Osher, 1985; Steger and Warming, 1981; Harten 
et al 1983) can be utilized. In the present work, the method of Roe (1981) was used in 
differencing the convective terms by an upwind method that is biased by the signs of 
the eigenvalues of the local flux Jacobian. This is accomplished by casting the governing 
equations in their characteristic form and then forming the differencing stencil such that 
it accounts for the direction of wave propagation. In the current formulation, the set 
of numerical equations are solved using a nonfactored line relaxation scheme similar to 
that employed by MacCormack (1985). This implicit scheme described in the next section 
makes use of large amounts of processor memory (approximately 180 times the number of 
grid points) to be coded in an efficient manner. 


3.6 Time- Accurate Algorithm Using Upwind Differencing 

Even though the current method is explained for time-accurate calculations, the 
same method can be used for steady state problems. 

3.6.1 Upwind Differencing 
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The upwind scheme for the convective flux derivatives is derived from 1-D theory, 
and then is applied to each of the coordinate directions separately. Flux- difference splitting 
is used here to structure the differencing stencil based on the sign of the eigenvalues of the 
convective flux Jacobian. The scheme presented here was originally derived by Roe (1981) 
as an approximate Riemann solver for the compressible gas dynamics equations. 

The derivative of the convective flux in the t direction is approximated by 

'*-1/2] 


where Ei+i j-i is a numerical flux and t is the discrete spatial index for the £ direction. 
The numerical flux is given by 

E i+1/2 = i [£(D i+1 ) + E(D.) - &+,/,] (3.57) 

where the is a dissipation term. For <^+ 1/2 = 0 this represents a second-order 

central difference scheme. A first-order upwind scheme is given by 

W = (a 1/2 - AEr +1/J ) (3.58) 

where AE± is the flux difference across positive or negative traveling waves. The flux 
difference is computed as 

A Ef +1/2 = 4 ± (5)AE i+1/2 (3.59) 

where the A operator is given by 

ADi + i/2 = Di+i — Di (3.60) 


dE [E{+ 1/2 — 

d( At 


(3.56) 


The plus (minus) Jacobian matrix has only positive or negative eigenvalues and is computed 
from 


A ± = XjAf X- 1 
A± = i(A I± |A,|) 


(3.61) 


where the subscript 1 denotes matrices corresponding to the ^-direction flux. The matrices 
X\ and X^ 1 are the right and left eigenvectors of the Jacobian matrix of the flux vector 
and Ai is a diagonal matrix consisting of its eigenvalues. All matrices appearing in the 
upwind dissipation term must be evaluated at a half point (denoted by i+l /2). To do this 
a special averaging of the dependent variables at neighboring points must be performed. 


33 



The Roe properties, which are necessary for a conservative scheme, are satisfied if the 
following averaging procedure is employed 

D = ±(D i+1 +Di) (3.62) 

A scheme of arbitrary order may be derived using these flux differences. Imple- 
mentation of higher-order-accurate approximations in an explicit scheme does not require 
significantly more computational time if the flux differences A E± axe all computed at once 
for a single line. A third-order upwind flux is then defined by 

<t>i+i/2 = -3 (& E t- j /2 - AEf +1/2 + A E7 +1/2 - A # i+3/2 ) (3.63) 

The primary problem with using schemes of accuracy greater than third-order oc- 
curs at the boundaries. Special treatment is needed, requiring a reduction of order or a 
much more complicated scheme. Therefore, when going to a higher-order-accurate scheme, 
compactness is desirable. Such a scheme was derived by Rai (1987) using a fifth-order- 
accurate, upwind-biased stencil. A fifth-order, fully upwind difference would require 11 
points, but this upwind-biased scheme requires only seven pointy. It is given by 

&+1/2 = — 3 q[ “ 2A Ef_ 3/2 + 11 A Ef_ 1/2 - 6AEf +1/2 - 3AE+ 3/2 ^ ^ 

+2A E~ +5/2 - 11AE~ +3 j 2 + 6A E7 +1/2 + 3A E7_ 1/2 ] 

Next to the boundary, near second-order accuracy can be maintained by the third- 
and fifth-order schemes by using the following 

. &+1/2 = e (AE+ +1/2 - AE7 + 1/2 ) (3.65) 

For e = 0, this flux leads to a second-order central difference. For e = 1, this is the same 
as the first-order dissipation term given by equation (3.58). By including a nonzero e, 
dissipation is added to the second-order, central-difference scheme to help suppress any 
oscillations. A value of e — 0.01 is commonly used for many applications. 


To form the delta fluxes used in this scheme, the eigensystem of the convective flux 
Jacobian is needed. For the current formulation, a generalized flux vector is given by 
equation (3.21) and the Jacobian matrix Ai = of the flux vector is given by equation 
(3.26). The normalized metrics Me redefined as 


k x 


ldtU 
J dx 


,* = 1,2,3 
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k - i - 1 2 3 
J dy' ' ' 

k - » - 1 2 3 

*‘ = 7ll ,< = 1,2,3 

As explained in conjunction with the diagonal algorithm in the previous section, a similarity 
transform for the Jacobian matrix is introduced here 


Ai = XiAiXr 1 


where A{ and c are defined by equations (3.39) and (3.40) respectively. The matrix of the 
right eigenvectors is given by 


where 


Xi = 


x k = 


Vk = 


Zk = 


Xfcfc = 


J Ikk = 


Zkk = 


0 

0 

/3(c - k t f 2) 

— 0 (c + kt/ 2 ) ' 


* 

u\$ + /3k x 

tiA4 + j3k x 

y* 

Vkh 

v\z + 0k y 

V^4 + /3ky 

Zk 

Zkk 

u>A 3 + (3k z 

10 X 4 + /3k z . 


dx 

df»+i 

dy 

dz 

d£*+i 

dx 

d£i + 2 
dy 

d£i+z 

dz 

dU+2 


&+i = Vi C) or £ for i = 1,2, and 3, respectively 


(3.66) 


(3.67) 


£*+2 = C>£> °r V for * = 1,2, and 3, respectively 
and its inverse can be obtained similarly. It should be noted that this transformation is 
nonsingular in any combination of metrics. 


3.6.2 Implicit Scheme 

This section describes the way in which equations (3.54) and (3.55) are numerically 
represented and solved. The first consideration is the formation of the Jacobian matrix 
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of the residual vector R required for the implicit side of the equation. Applying the 
difference formula given in equation (3.56) to the convective flux vectors and applying a 
second-order, central- difference formula to the viscous terms, the residual at a discrete 
point (xij k ,yijk,Zijk) is given by 


Rijh = 


Ei+l/2,j,k ~ Ei-l/2,j,k E itj+1 /2,k - Ejj- 1/2, fe fc+1/2 - 1/2 


AC 


+ 


At} 


+ 


AC 


_ (Ey)i+ l,j,k ~ (Ey)i- ~ {Ey)i,j- l,fc _ 

2AC 2A77 2A77 

(3.68) 

The generalized coordinates are chosen so that AC, Atj, and AC axe equal to one. To limit 
the bandwidth of the implicit system of equations, the Jacobian of the residual vector will 
be formed by considering only first-order contributions to the upwind numerical fluxes as 
well as the second-order differencing of the viscous terms. Hence, the only portion of the 
residual vector that is actually linearized is the following. 


Rijk — Ei-l t j t k + ... 

~^ E t+l/2,j,k + ^ E i+l/2,j,k + AE i-l/2,j,k ~ A E i-i/ 2 ,j,k 

-(E v )i+i/2,j,k + (Ev)i-l/2,j,k - — 


(3.69) 


The exact Jacobian of this residual vector will result in a banded matrix of the form 


dR 

dD 


= B 


dRa 


■ijk 


9Dij,k-i 


,0,...,0, 


dRjjk 
dDij-i'k 


-, 0, ..., 0, 


dR ijk dRjjk dRjjk 

dDi-ij t k ’ dD itj ,k ' dD i+ltj> k ’ 


0, 


0, 


dRjjk 

dDjj+i'k 


, 0 , 


0, 


dRjjk 

dDjj t k+i . 


(3.70) 


These exact Jacobians can be very costly to form. Therefore, approximate Jacobians of 
the flux differences sis derived and analyzed by both Yee (1986) and Barth (1987) are used. 
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These are given as follows 


dRjjk 

dD itjtk -i 

d Rjjk 

dD iJ)k 


W 2^ Cj t j,k-1 ^t,j,k-l/2 + ^i,j,k-\/l) (T3 1/2 

K 2 (^t+l/2,j,k + ^t-1/2, j,k ~ ^H-l/2,j,fc “ ^i-1/2 ,j,k 

+ B i,j+l/2,k + B t,j-l/2,k ~ B i,j+l/2,k ~ B i,j-l/2,k 
+ C i,j,k+ 1/2 + C t,j,k-\/2 ~ C i,j,k+ 1/2 ~ C i,j, k-l/ 2 ) 

+ (7l)*+l/2,j,fc + (72)*,j+l/2,fc + (73)i,i,fc+l/2 
+ (7l)i— 1/2, j, As + (72 )i,j — l/2,fc + (73)t,j,fc-l/2 


dRjjk 

dDj+\,j,k 


W 2^ i + 1 'i' 1 -^i+l/2,j,fc + ^*+l/2,i,Jfe) (7l )i+l/2,j,fc 


(3.71) 


where A = Ai,B = A 2 , and (7 = A 3 as given by equation (3.26), and where only the 
orthogonal mesh terms are retained for the implicit viscous terms. 

The matrix equation is solved using a line-relaxation method. First, the entire 
numerical matrix equation is formed from values at the nth time level. A significant 
amount of memory is required to store all of these values, but it is justified by the savings 
in computing time. At this point the numerical equation is stored as the following banded 
matrix 


B[T, 0, ..., 0, U, 0, ..., 0, X, Y, Z , 0, ..., 0, V, 0, ..., 0, W]AD = R 


(3.72) 


where A D = D n+1 — D n and T y U,V,W,X,Y, and Z are vectors of 4 x 4 blocks which 
lie on the diagonals of the banded matrix, with the Y vector on the main diagonal. This 
matrix equation is approximately solved using an iterative approach. One family of lines 
is used as the sweep direction. Using, for example, the £ family, a tri diagonal matrix is 
formed by multiplying the elements outside the tridiagonal band by the current AD and 
shifting them over to the right-hand side. This system can be represented by the following 


B[X, 7, Z]AD l+1 =R- TAD 1 ^, 

~ VAD\ tj+lth 


~ UAD l itj _ ltk 
~ WAD \ tj ' k+1 


where l is an iteration index. The number of iterations is generally very small. This equa- 
tion can be solved most efficiently by first performing and storing the LU decomposition 
of the tridiagonal matrix before the iteration is begun. Then for each iteration, the right- 
hand side is formed using the latest known AD (which is set to zero for the first iteration), 
and the entire system is backsolved. The LU decomposition can be entirely vectorized, 
but the back solution is inherently recursive and cannot be vectorized. 
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3.6.3 Boundary Conditions for Upwind Scheme 


Implicit boundary conditions are used at all of the boundaries, which enables the 
use of large time steps. At a viscous no-slip surface, the velocity is specified to be zero, and 
the pressure at the boundary is obtained by specifying that the pressure gradient normal 
to the wall be zero. The boundary conditions used for the inflow and outflow regions are 
based on the method of characteristics. The formulation of these boundary conditions 
is similar to that given by Merkle arid Tsai (1986), but the implementation is slightly 
different. The scheme is derived here for a constant £ boundary, with similar results for a 
constant 77 or a constant ( boundary. The finite-speed waves which arise with the use of 
artificial compressibility are governed by the following 

dD dE _ dEdD , dD _ YAY _i dD 
dr ~ di~ dD di~ di~ di 


Multiplying by X 1 gives 




(3.73) 


If the X~ l matrix is moved inside the spatial- and time-derivative, the result is a 
system of independent scalar equations, each having the form of a wave equation. These 
equations axe termed the characteristic equations. The sign of the eigenvalues in the A 
matrix determines the direction of travel of each of the waves. For a positive or negative 
eigenvalue, there corresponds a wave that propagates information in the positive or negative 
£ direction. The number of positive or negative eigenvalues determines the number of 
characteristic equations propagating information from the interior of the computational 
domain to the boundary. Thus, at the boundary, the characteristics equations which bring 
information from the interior will be chosen as part of the boundary conditions. The rest 
of the information should come from outside the computational domain, which leaves some 
variables to be specified. For inflow boundaries, two different sets of specified variables 
have been used successfully. One set consists of the total pressure and the cross-flow 
velocity. This set is useful for problems in which the inflow velocity profile is not known. 
For outflow boundaries, static pressures have been specified for computations presented 
later in this note. 
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IV. A FRACTIONAL STEP METHOD IN GENERALIZED COORDINATES 


In the previous chapter, a pseudocompressibility method was described where a 
fictitious time derivative of pressure is added to the continuity equation so that the modified 
set of equations can be solved implicitly by marching in time. When a steady solution 
is reached the original equations axe recovered. Time accuracy could be achieved by 
employing an iterative technique in a pseudo time level. An alternative approach for time- 
dependent computations is the pressure iteration method as outlined in chapter II. Among 
several choices in this category, the fractional step method allows flexibility in combining 
various operator splitting techniques. In light of successful computations in Cartesian 
coordinates using its numerous variants, we were interested in developing a general 3-D 
solver based on a fractional step approach. In this chapter, some details of the method 
developed by Rosenfeld et al. (1988, 1989) will be discussed. 

4.1 Overview 

Several key features must be carefully handled in developing this flow solver. Accu- 
racy in an arbitrary domain is related to the method of discretizing the governing equations. 
Geometric quantities, for example, have to be discretized consistently with the flow differ- 
encing. For computational efficiency, an efficient Poisson solver for pressure is critical in 
this approach. To accomplish this, variable definitions and mesh arrangement are chosen to 
develop an efficient Poisson solver in curvilinear coordinates. A ZEBRA scheme with four 
color ordering has been devised for the present Poisson solver, in which a multigrid accel- 
eration procedure can be readily incorporated in the future. In this respect, the governing 
equations are discretized by a finite-volume scheme on a staggered grid. As for the choice 
of dependent variables, the volume fluxes across each face of the computational cells are 
used. Thus the discretized mass conservation equation can be satisfied consistently with 
the flux balancing scheme explained later. In solving the momentum equations, a sig- 
nificant part of the viscous terms axe solved implicitly by an approximate factorization 
scheme to minimize time step limitation resulting from use of a viscous grid. A staggered 
grid has favorable properties in cartesiati coordinates, such as coupling odd-even points. 
However, it is debatable whether a staggered grid has clear advantages over a regular grid 
in generalized curvilinear coordinates. In the present procedure, a staggered arrangement 
is chosen to simplify the pressure boundary condition in devising a Poisson solver. 
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4.2 Formulation 


The equations governing the flow of isothermal, constant density, incompressible, 
viscous fluids in a time- dependent control volume with the face area vector S(<) and volume 
V(t) are written in integral form for the conservation of mass 

dV r 

— + j> s dS-(u-v) = 0 (4.1) 

and for the conservation of momentum 

If 

dt Jv 

where t is the time, u is the velocity vector, v is the surface element velocity resulting 
from the motion of the grid, and dS is a surface area vector. The tensor T is given by 

T = — (u — v)u — pi + v (Vu + (Vu) T ) (4.3) 

where Vu is the gradient of u while (*) T is the transpose operator. 


u dV 


=1 


dS-T 


(4.2) 


The only differences between the fixed and the moving-grid equations are the terms 
that include the surface element velocity v and the time dependence of the cell geometry 
(volume and face area). The volume conservation of each time- varying cell requires 


8V 

dt 



d S • v = 0 


(4.4) 


where the term dS • v represents the volume swept out by the face S over the time increment . 
dt. Thus, the mass conservation equation has exactly the same form as for fixed grids 



dS • u = 0 


(4.5) 


The usual practice is to transform equations (4.2) and (4.5) into a differential form. In 
the present work, the integral formulation is maintained for convenience in deriving the 
finite- volume scheme for arbitrarily moving geometries. 
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4.3 Discretization 


In the present integral formulation, spatial derivative terms are discretized in a 
different way from the differential equation approach as used in the previous chapter. In 
this section, details of the spatial discretization method are described. 


4.3.1 Geometric Quantities 

A general nonorthogonal coordinate system (£,*7, C) is defined by 


(4.6) 


The center of each primary cell is designated by the indices (i,j,k). The area of the face 
I of a primary cell (figure 4.1) is given by the vector quantity 

_7 dr dr 

S ~ 9(1 + 1) X 9(7+2) ( ’ 7) 

where the computational coordinates I = £, rj or £ are in cyclic order and x is the cross 
product operator. The vector S l has the magnitude of the face area and a direction normal 
to it and is the contravariant base vector VI scaled by the inverse of the Jacobian 1/ J, i.e., 
S l = j VI. As pointed out by Vinokur (1986), am accurate discretization should satisfy 
certain geometric identities. The condition that a cell is closed should be satisfied exactly 
in the discrete form 

£ S' = ° (4.8) 

l 

where the summation is over all the faces of the computational cell. Equation (4.8) cam be 
satisfied if S* is approximated from equation (4.7) by a proper approximation of |y. For 
example, is computed at the point (* + 1 , j, k) by using the second order approximation: 

dr 1 

^dr^ i+ \ = ~ T j-b h ~h +r i+£.*+£ - r j-£.H-§)*+i 

( dr v 1 , N 

"*' r i+5.*+i ~ r J+f. fe -^*+5 

The volumes of all discrete computational cells will sum up to the total volume at a 
given time if the volume of each computational cell is computed by dividing the cell into 
three pyramids having in common the main diagonal and one vertex of the cell, resulting 
in the following 


V = 5 ( s f-i + S J-J + S Lj) • (r i+iJ+M+ } -K-JJ-i.*-}) 


(4.9) 
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In the present method, the volume conservation equation (4.4) is satisfied discretely by 
interpreting the term dS • v in equation (4.4) as the rate of the volume swept by the face 
dS. For example the volume swept by the face 1 can be computed by a formula similar 
to equation (4.9). 


+ SS U + (4 ' 10) 

where the time level is given by n. The quantities 6S7, and axe the areas swept 

by the motion of the face S^,_ i - see the shaded areas in figure 4.1b. The area 6 S : t is 

J » j u j — j 

computed from 

6S U - 5 - r ‘+i> x (r 2» - (4 ' n) 


and SSl_ ^ can be computed in a similar way. The volume of the cell at the time level 
(n + 1) is computed from equation (4.4). 


yn+i _ V n + ^(^) n +i (4.12) 

i 

Note that 6V l /At is the volume flux resulting from the motion of the coordinate system 
and has a meaning similar to the volume-flux U l . An accurate computation of the volume 
of each cell at the time level (n + 1) is important for the accurate representation of the 
momentum equation. 


In the present finite volume formulation, no coordinate derivatives appear directly 
in the discrete equations, as in the case of finite- difference formulas. Instead, quantities 
with clear geometric meaning such as the volume and the face areas of the computational 
cells are used. The discrete approximation of these quantities satisfies the geometric con- 
servation laws. A principal difference between the finite-volume and the finite-difference 
approach to moving grids is in the interpretation of the quantity 8V l j At. In the finite- 
volume method, it is treated as a geometric quantity which expresses the rate of displace- 
ment of a cell face, whereas in the finite-difference method, the grid velocity is combined 
with the fluid velocity to define a “relative flow velocity” (see Vinokur, 1986). 


4.3.2 Mass Conservation Equation 

The discretization of the mass conservation equation (4.5) over the faces of the 
primary computational cells yields 

(S* • u) i+ i - (S« • u)i_i + (S” • u) i+ i - (S” ■ u)^! + (S< • u) fc+ i - (S< • u) fc _i = 0 (4.13) 
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Note that the default indices and the time-level (to + 1) axe usually omitted for 

simplicity. Each term on the left-hand side of equation (4.13) approximates the volume- 
flux over a face of the primary cell. A simple discretized mass conservation equation can 
be obtained by using the following variables 

TJ* = S* • u 

U v = S’ • u (4.14) 

U< = S c • u 

as the unknowns instead of the Cartesian velocity components. The quantities U^,U V , 
and are the volume fluxes over the £, rj, and ( faces of a primary cell, respectively. In 
tensor algebra nomenclature, these are the contravariant components of the velocity vector 
scaled by the inverse of the Jacobian (1 /J). With this choice of the dependent variables, 
the continuity equation takes a form identical to the Cartesian case 

< j - U f-i +U] +i -U]_ i + Ui +i - Ui_ J = XW7* = 0 (4.15) 

As explained in chapter II, it is crucial to satisfy the discrete mass conservation 
equation. Therefore, the simple form of equation (4.15) suggests that the volume fluxes 
can be chosen as the ‘natural’ dependent variables for fractional step methods. On the 
other hand, this choice complicates the discretization of the momentum equations. In the 
present work, volume fluxes and the pressure are chosen as the dependent variables. 

4.3.3 Momentum Conservation Equation 

Spatial discretization of the momentum-conservation-law equation (4.2) for a com- 
putational cell with volume V yields 

i(Vu) = ^S'-f = F (4.16) 

Following the development of the scheme given by Rosenfeld et al. (1988), a general 
two-step temporal discretization of equation (4.16) for a constant time-step A t is 

(l+e)(Vu) n+1 — (l+2e)(Vu) n -f-e(Vu) n-1 = A<[0F n+1 +(l-^)F n +<^(F n -F n - 1 )] (4.17) 

where to is the time level and e, 0, and <f> are arbitrary parameters. A truncation error 
analysis shows that second-order accuracy in time is obtained if 

^ + e = 6 + <f> (4-18) 
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In order to replace u by the new dependent variables U l , the corresponding area vectors 
are dotted with the momentum equations and the integral momentum equation is ap- 
plied on different computational cells for each unknown U l . Each cell has the dimensions 
At x Arj x A£, but the centers Eire located at ( i + |,y, &), (i,j + |,fc) and |) for 

the U^ y U v , and momentum equations, respectively. 

The derivation of the ^-momentum equation will be described next. The other two 
momentum equations can be obtained by cyclic permutation. Using the identity 


u = S tU* +S V U V + S C U C = S m U m 

the dot product of equation (4.17) and results in 

(1 + e)(V U*) n+1 - (1 + 2e)(V U l ) n (S< • S?) + e(V U l ) k -\ S* • Sf _1 ) 
= At[$S* ■ F” +1 + (1 - 6) S ( • F n + <f>S ( • (F n - F n-1 )] 


(4.19) 


(4.20) 


where S m is the inverse base of S*. Note that the n+1 is the default time-level and therefore 
S* is computed at n + 1 . To minimize the temporal truncation error, the coefficients in 
the above equation are chosen as 


t= \' ^ = 


0 = 1 


Then equation (4.20) becomes 

3(V P<)" +I - i(V U‘) n ( S< ■ Sf) + ( V J7')" _I (S f • S" _I ) = 2A<I" +1 


(4.21) 


where 


L = S* • F 


The difference between a fixed-grid and a moving-grid case is in the computation of the 
convection terms which should include the motion of the grid. For example, the convection 
flux of the ^-momentum equation on the £-face center (»,j,n) is given by 




The difference equations are second-order accurate in time if 6V* would not lag in time 
by ^ over the volume-flux terms U 1 (see eq. (4.10)). The resulting discrete equations are 
conservative in any moving coordinate system and are spatially second-order accurate. For 
high Reynolds number flows, fourth-order dissipation is added to eliminate high frequency 
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components of the solution. The dissipation terms are interpreted in terms of fluxes and 
therefore the conservation properties of the equations are maintained. 

The above time-marching scheme requires the computation of the operator L (. 

£ < = s f+i • E s ' • f < 4 - 22 ) 

where the summation is over all the faces of a computational cell. It is to be noted that 

U l = S' • u = S' • S m U m (4.23a) 


and the invariance of the velocity vector requires: 

S' • S m = S lr 


(4.236) 


where $( m is the Kronecker delta, and S 1 is the inverse base to S m and has the differential 
dr 

analogue S m = J-z — . In terms of tensor algebra, S m is the covariant base vector scaled 
(yrrt 

by the Jacobian J while S 1 is the contravariant base vector scaled by 1/J. A uniform 
velocity field can be numerically preserved if the base S m is computed at each point from 
the relation 


Sm — 


gm-f-1 ^ gm-f2 


(4.24) 


S m • (S m+1 x S m + 2 ) 

which satisfies (4.23b) identically. The variable m is the cyclic permutation of (£, *7, C)- 

The product S' • T should be computed for each face of each momentum equation, 
see equation (4.22). For example, the ( face-center for the U* momentum cell is located 
at (t,j,n) (see fig. 2). The flux over this face is computed from 

(S< ■ T) iJtn = (~{U< - ^)U‘ S, - S (p + S< • v(Vo + (Vu) T )) (4.25) 


The conservative form of the velocity vector gradient is 

Vu 


_ $s dSu 


(4.26) 


Applying equation (4.26) for the computation of Vu» (J|n yields 



£-4 u ‘ + j 



(4.27) 
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where only those indices different from ( i,j,n ) are given. 

The rj or ( face-centers are located at (*+ |, y — j,n) and (i+ 5 , j, n— |), respectively. 
The fluxes over these faces are computed in a similar way. The convection and diffusion 
fluxes in equation (4.25) may be approximated in various ways. In the present work, all 
of the unknowns at the point (t, j, n) are computed by simple averaging and therefore the 
scheme is equivalent to the results computed by second order central differencing. 


4.4 Solution Procedure 


The dependent variables are the pressure, defined at the center of the primary 
cells, and the volume fluxes defined on the faces of the primary cells. This selection 
is equivalent to a finite- difference formulation over a staggered grid with the choice of 
scaled contravariant velocity components as the unknowns. The mass and momentum 
conservation equations (4.5) and (4.20) are solved by a fractional step method. 

First the momentum equations are solved for an approximate U l . This is done by 
computing AU l defined by 

A U l = ( U l ) n+1 - ( U l ) n 

This first step of the velocity correction is obtained from equation (4.20) by evaluating the 
pressure gradient term at time step (n). As explained in chapter II, various time advancing 
schemes can be implemented here. One possibility is an Adams- B as hforth type scheme 
used by Rosenfeld et al. (1988) which can be written as follows 

(v; + i I - A t.Di) A U l = A t Qi? - - i^(AP n " 1 ) - D^AU^- 1 )^ (4.28) 

where Di represents the viscous terms and Ri the pressure gradient terms. Here, only 
orthogonal viscous terms are included on the implicit side. Later linearized convection 
terms were included implicitly (Rosenfeld and Kwak, 1989), thus a larger time step could 
be taken. Discretization of each term in the above equation is given in the previous section. 


In the second step, the velocity is updated such that the continuity equation is 
satisfied at the next time level. This is achieved by a single Poisson-type equation as 
follows: 



(4.29) 
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As explained in chapter II, it is to be noted that the Laplacian operator is written as a 
discrete divergence of a discrete gradient operator. This is an important step to satisfy 
the mass conservation equation in discretized space. Then the variables at the new time 
level n + 1 are computed by 

(i U l ) n+1 = (U l ) n + AU l 

pn+l = pn + ^ 



It can be shown that U l k+1 is the exact discrete solution of the momentum and mass 
conservation equations, but A P — O(Ai) and therefore P n+1 is not the exact discrete 
solution. In a Cartesian coordinate, Kim and Moin (1985) used the exact relationship 


P" + ‘ = *- 


where V 2 is the Laplacian operator. However, such a simple relationship cannot be found 
in generalized coordinate systems with the present splitting of the diffusion terms. Braza 
et al. (1986) have not found significant differences between solutions which used the 
approximation P = <j> and solutions which computed the discrete pressure exactly, although 
the approximate solution may be less stable. As the difference between (f> and AP is 
proportional to v, the direct substitution AP = <f> is reasonable for high Reynolds number 
flows. However, this approximation may degrade the formal second order accuracy in time 
of the present method. This issue requires additional study. 


The first step is a consistent approximation of the momentum equations, e.g. as 
(At, Al) — + 0, Z7j[ +1 — ► U l , where U l is the solution of the continuous problem. Therefore, 

the physical boundary conditions may be specified at the (n + 1) time level. In some 
fractional step methods, the steps are not consistent and special boundary conditions 
should be devised for the intermediate steps. 

One drawback of the present method is the increased memory required for storing 
variables at three time levels. Also, the effect of the many geometric approximations on 
the overall accuracy is yet to be fully evaluated. Even though the current procedure has 
been validated by many test cases including the moving domain problem, the relative 
advantages or disadvantages to pseudocompressibility method will have to be determined. 
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V. DEVELOPMENT AND VALIDATION OF FLOW SOLVERS 


In the previous two chapters, solution methods based on two major primitive vari- 
able approaches, namely, pseudocompressibility and a fractional step approach, have been 
described for 3-D applications in generalized coordinates. Several computer codes have 
been developed according to these procedures. Once these codes are written, the accuracy, 
convergence speed, computational efficiency, and robustness have to be carefully evaluated. 

The code development effort was motivated by a demand for efficient flow solvers 
in severed major application areas, namely, simulation of the flow through an impeller in 
the SSME power head and other complex internal flows, simulation of a full aerodynamic 
configuration at low speed involving a wing-body type juncture flow, and time-dependent 
flow analysis of an artificial heart. Validation of the above flow solvers was done by 
choosing several basic configurations relevant to the local flow encountered in the above 
applications. To remove uncertainties caused by turbulence modeling, validation cases 
were chosen in the laminar range of Reynolds numbers first followed by turbulent cases. In 
this chapter, a summary of these codes and several representative validation computations 
will be presented. 


5.1 Code Development 

The following is a list of these incompressible Navier-Stokes codes. Each code is 
given a name and the results presented here will be identified by these names. 

5.1.1 INS3D Code 

A computer code hasbeen developed based on the steady-state algorithm described in sec- 
tion 3.4. A steady formulation of the pseudocompressibility approach is solved using an 
approximate factorization scheme. This takes advantage of the recent advances in state- 
of-the-art CFD made in conjunction with compressible flow computations. The spatial 
discretization utilizes second-order central differencing with additional numerical dissipa- 
tion terms. The code has been validated and numerous 3-D problems have been solved 
using this code. 
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5.1.2 INS3D-UP Code 


To obtain time-accurate solutions using the pseudocompressibility formulation, it is neces- 
sary to satisfy the continuity equation at each time step by subiteration in pseudotime. In 
order to use a large time step in the pseudotime iteration, an upwind differencing scheme 
based on flux-difference splitting is used combined with an implicit line relaxation scheme. 
This removes the factorization error and the need for numerical dissipation terms. The 
code has been validated. Major unsteady flow applications include simulation of the flow 
through an artificial heart. 

5.1.3 INS3D-LU Code 

This code is also based on the pseudocompressibility formulation. However, a finite vol- 
ume scheme in conjunction with either central or upwind differencing is used for spacial 
discretization. An LU-SGS implicit algorithm is employed for temporal discretization. 
This code has an option to utilize a rotating coordinate system so that rotor-steady-state 
solutions (steady in relative frame of reference) can be computed. Therefore, the flowfield 
of a propeller in an infinite medium or in a ducted enviroment can be solved using this 
code. The code is completely vectorized; validation is in progress in conjunction with the 
Space Shuttle main engine flow analysis. 

5.1.4 INS3D-FS Code 

A generalized flow solver based on a fractional- step method has been developed for time- 
dependent computations of the incompressible Navier-Stokes equations. The governing 
equations are discretized conservatively using finite-volume approach on a staggered grid. 
Here, the discretized equations axe advanced in time by decoupling the solution of the 
momentum equation from that of the continuity equation. This procedure, combined 
with accurate and consistent approximations of the geometric quantities, satisfies the dis- 
cretized mass conservation equation exactly in a discrete sense. A four-color ZEBRA 
scheme is devised for solving the Poisson equation for the pressure correction. Several 
sets of computational results have been compared with experiments and other numerical 
solutions. 


5.2 Channel Flow 

The physical interpretation of the pseudocompressibility was given in chapter III. It 
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is therefore interesting to test the validity of the estimate on the permissible range of the 
pseudocompressibility parameter, (3. To do this, a validation calculation is performed over 
a range of /? using a simple test problem. This study is done using the INS3D code. The 
channel flow is perhaps the simplest internal-flow test problem, where the pressure wave 
propagates between the in- and outflow boundaries while the viscous effect spreads inward 
from two walls. The coordinate system of a 2-D rectangular channel with a width of 1 
and length of 15 is shown in figure 5.1, where velocity vectors for a converged solution are 
also shown. To obtain fully developed velocity profiles within a reasonable channel length, 
a partially developed boundary layer profile can be imposed at the inflow boundary. In 
this experiment, a uniform inlet flow is prescribed, and the Reynolds number based on 
the duct width and the average velocity is 1000, and the time step At = 0.1. Then the 
recommended range of (3 is estimated to be 

0.12 <£<10 

To illustrate the pressure wave propagation phenomena and its effect on the convergence 
property, the channel is impulsively started. Here, five different values of f3 (0.1, 1, 5, 10 
and 50 ) were chosen such that two cases axe outside the range while three valuesare kept 
within the range. 

Pressure contours at three different time levels are shown in figure 5.2, for /? = 5 
which is within the recommended range. The expansion wave from the exit plane propa- 
gates upstream sufficiently fast to balance the spreading of the viscous effect. The solution 
converges nicely in this case. However, for /3 = 0.1 which is lower than the recommended 
lower bound of /3, the speed of the upstream propagating pressure wave from the exit plane 
is very low. Therefore, the expansion wave is confined near the exit plane while the viscous 
effect spreads into the flow field from both upper and lower surfaces. Thus the viscous 
field is not properly balanced by the right pressure gradient which causes the spurious 
fluctuations to amplify, as shown at r = 2.5, and eventually the computation blows up at 
r = 4.0. 

The convergence history is shown in figure 5.3(a). The log of the root-mean-square 
of the change in pressure and velocities (RMSDQ) is plotted against the computation time 
r. It can be seen that calculations for f3=Q.l and 50 become unstable within 50 steps and 
start to diverge, whereas other cases converge to a stable solution. The effect of j3 values 
on the incompressibility of the fluid is shown in figure 5.3(b) in the form of the log of 
the root-mean-square of the divergence of the velocity field (RMSDIV) plotted against the 
pseudo time r. 
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5.3 Flow Over A Backward-Facing Step 


The flow over a 2-D backward-facing step is computed using INS3D and INS3D- 
UP. Figure 5.4 shows the schematic of the problem, where the step height is equal to the 
inlet height and the Reynolds number is based on twice the step height. The upstream 
boundary is located at the step and a fully developed channel flow is imposed at the inlet. 
This problem is very challenging computationally as it involves a primary and a secondary 
separation bubble. The size and the location of these separation zones are very sensitive 
to the pressure gradient, thereby providing a good viscous flow validation case. 


Results obtained using INS3D have been reported previously by Rogers et al. (1985). 
The separation lengths were found to be sensitive to the magnitude of the numerical dissi- 
pation coefficient at higher Reynolds numbers. This in essence is equivalent to changing the 
effective Reynolds number of the flow as the Reynolds number approaches 1000. Therefore, 
it is quite desirable to remove the dissipation model dependence on the solution at a high 
laminar range of Reynolds numbers. In light of the similarity of the pseudocompressibility 
formulation to that of the compressible flows, the study done by Yoon and Kwak (1988) 
may be applied to achieve this. When the flow becomes turbulent, the Reynolds number 
based on turbulent eddy viscosity is in the order of a few hundred, and thus the present 
dissipation model may perform adequately. 

More extensive validation has been done by using the INS3D-UP code by Rogers 
(1988). The computed separation and reattachment locations are then compared to exper- 
imental values given by Armaly et al. (1983) in figure 5.5. For the primary reattachment 
length, xl, good agreement is observed between the experiments and the computation 
until the secondary separation appears at a Reynolds number of about 400. At a higher 
Reynolds number, the primary separation length, xl, and the secondary separation point, 
x2, deviates from the experimental data. Similar results were reported by Kim and Moin 
(1985) using a fractional step method. Armarly et al. reported that 3-D flow was observed 
near the step when the Reynolds number is greater than 400. For more complete validation, 
the numerical simulation needs to include the entire experimental configuration including 
the possible 3-D effect. In figure 5.6, the convergence history is shown. Using the grid 
of 100 points in the streamwise direction and 53 points in the crossflow direction, a good 
converged solution is obtained within 100 iterations and under 11.5 seconds of computing 
time on the Cray 2. 
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5.4 Flow Over A Circular Cylinder 


To represent an external flow, the flow past a circular cylinder was computed. For 
external flows in which the computational domain extends a large distance from the body, 
the pressure waves originating from the body surface propagate into the farfield. Therefore, 
to obtain the near-field solution only, the distance traveled by the waves and the spreading 
of the vorticity can be considered approximately the same in magnitude. The range of 
0 in INS3D can then be estimated based on this reasoning. For example, at Re=40, if 
the viscous region is taken to be approximately two diameters away from the body, the 
following range for 0 using At = 0.1 can be estimated as 

0.1 « 0 < 10 . 

Physically, for external flows, the pressure wave can quickly travel a short distance to 
balance the viscous region close to the body. Therefore, the magnitude of 0 is less restrictive 
than for internal flow cases such as the channel flow described above, which is consistent 
with our past experience. Good agreement with literature data was found for steady state 
computation at Re=40 (Kwak et al. 1986, INS3D). Far more extensive validation was done 
at Re=5, 10, 20, and 40 by Rogers and Kwak (1988) using INS3D-UP and showed good 
agreement with experiments and other computations (not illustrated here). 

To capture the near-field detail of transient flow, impulsively started circular cylin- 
ders at Re=40 and 200 is computed using INS3D-FS code (see Rosenfeld et al., 1988, 
1989). In figure 5.7, time evolution of the separation length is compared with experiments 
by Coutanceau and Bouard (1977) and other computations by Collins and Dennis (1973). 
The comparison is very good. 

As the Reynolds number increases above 40, nonsymmetric wake develops and pe- 
riodic vortex shedding sets in (for a comprehensive review, see Morkovin, 1964). Both 
INS3D-UP and INS3D-FS codes are validated using this simple yet challenging problem. 
These computations are compared with other numerical and experimental results. The 
calculation was performed using an 0-type grid clustered near the body. To obtain time- 
accurate solutions from INS3D-UP, the subiterations were carried out at each physical time 
step. Starting impulsively from rest, over 20 subiterations were required during the tran- 
sient phase. Non-symmetric wake develops spontaneously followed by shedding of vortices 
without introducing any artificial disturbance, probably due to biasing in the upwinding 
scheme which in turn may have introduced enough disturbance into the flow field. When a 
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central differencing scheme is used as in INS3D-FS, it is necessary to introduce asymmet- 
ric disturbance to initiate the shedding. Figure 5.8 shows a Strouhal number plot versus 
Reynolds number from computed results compared to experiments. Good agreement is 
observed. The Strouhal number is perhaps relatively easy to predict. However, the time- 
accuracy of the code itself has been validated using other problems with exact solutions 
(see Rogers and Kwak, 1988a). 

The staggered pattern of the vortex shedding known as Karman’s vortex street 
has been a subject of many flow visualization studies. For the purpose of comparison, 
particle traces are generated from the time-dependent solution of flow around a circular 
cylinder at a Reynolds number of 105 using the INS3D-UP (see Rogers and Kwak, 1988). 
Figure 5.9(a) shows this computed vortex street. Shown in figure 5.9(b) is an experimental 
photograph of the same conditions taken by Taneda and reproduced from Van Dyke (1982). 
The streaklines in the experiment are shown by electrolytic precipitation in water. The 
vortex structure is seen to be very similar between the two. The experimental picture is 
digitized and displayed on a workstation along with the computationally generated flow 
visualization picture. This example illustrates the potential of using postprocessing of the 
CFD results for studying fundamental fluid dynamic phenomena. 


5.5 Curved Duct of Square Cross-Section 

The flow through a square duct with a 90° bend offers a good test case for a full 
Navier-Stokes solver. This flow is rich in secondary flow phenomena both in the corner 
regions and through the curvature in a streamwise direction. This geometry was studied 
experimentally by Humphrey et al. (1977) and Taylor et al. (1981, 1982), and extensive 
laminar flow data are available. This particular geometry was used as a steady state test 
case for both INS3D and INS3D-UP code. The geometry is shown in figure 5.10; the 
Reynolds number of the flow is 790. The problem was nondimensionalized using the side 
of the square cross-section, H. 

For validating INS3D-UP (Rogers and Kwak, 1989), the inflow velocity was specified 
to be that of the fully developed, laminar, straight square duct (see White, 1974). The 
velocity is normalized by the average inflow velocity. Three different grids were used and 
the pseudocompressibility parameter 0 wets set to be unity throughout the computation. 
The computed results are compared to the experimental results of Humphrey et ad. (1977) 
as shown in figure 5.11. Overall, the comparison is quite satisfactory. The velocity profile 
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through the middle of the bend is very good. However, at the latter part of the bend, the 
computation shows a double inflection in the velocity profile, which is not evident in the 
experimental results. It is seen that the peak velocity is very well predicted. 

Extensive validation of the INS3D code was performed by McConnaughey et al. 
(1989). Utilizing the detailed inflow measurement by Taylor et al. (1981, 1982), they 
investigated various aspect of the flow solver including a grid refinement study and a 
sensitivity study of the numerical dissipation terms. Computed results on a 90° bend 
and an S-bend are extensively compared in the same report. Figures 5.12 and 5.13 show 
detailed velocity contours which illustrate the quality of the numerical solutions compared 
with experiments. The flow field details compare very nicely and the overall comparison 
is excellent. 


5.6 Juncture Flow 

The flow around a cylinder-plate or a wing-body juncture produces an interesting 
viscous phenomena due to the interaction between the boundary layer from the plate 
and the cylinder. The 3-D separation of the boundary layer and subsequent formation 
of the so-called horseshoe vortex and its development is very challenging to analyze both 
experimentally and numerically. This type of flow can occur in many practical engineering 
problems. Flow around a wing-fuselage junction and an appendage-submarine body are 
simple examples of its kind, and the flow near the endwall of turbomachinery blades mi g ht 
be one of the most complicated juncture flow problems. One major motivation of studying 
this type of juncture flow is related to the flow analysis of the Space Shuttle main engine, 
which will be presented in detail in the next chapter. In the SSME, liquid oxygen posts are 
densely packed in the main injector region. Even though a single cylinder-plate flow is an 
extreme idealization of the flow in the actual SSME environment, validating the computer 
code in this simplified case is of considerable value in extending the code to realistic case 
where detailed experimental measurement is very difficult and scarce. 

5.6.1 Cylinder-Flat Plate 

Most of the earlier studies on cylinder-flat plate juncture flow have been experimen- 
tal. Baker (1979) shows that laminar juncture flow is confined to a very limited region. A 
similar result has been obtained most recently by Thomas (1986). Eckerle and Langston 
(1986) reported a single primary vortex and saddle point contrary to multiple vortex sys- 
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terns observed earlier by other researchers. Interpretation of the phenomena also varies 
(Thomas, 1986 and Peak and Tobak, 1980). 

Computational simulation of these flows involves distinctively different features from 
those of external aerodynamics. For instance, the thickness of the viscous layer for these 
types of flows is of the same order as the characteristic flow-field dimension, while the 
viscous region tends to be confined in a thin layer near the body for external flows. Realistic 
juncture flows under an internal flow environment are likely to have a large amount of 
deflection as in the case of LOX post regions in the SSME. Recently, some numerical 
studies on this flow have been attempted. Kaul et al. (1985) reported a numerical study 
on a single cylinder- plate flow using INS3D code. Highlights of this and other studies 
were reported by Kwak et al. (1986). Independently, Kiehm et al (1986) reported a 
numerical study of flow around a single post in a channel. These computational results 
show qualitatively similar phenomena. Here, several representative results using INS3D 
are presented. 

In figure 5.14, the computational domain for a single post on a flat plate is illustrated. 
The upstream boundary layer thickness is varied by using partially and fully developed 
channel flow profiles. The convergence characteristics of the flow solver are shown in 
figure 5.15 by the history of RMSDQ, which denotes the root-mean-square value of the 
change per iteration in the pressure and velocities. The three curves in the figure show 
three variations of the INS3D code; namely a block tridiagonal, a diagonal version with 
second-order implicit smoothing terms, and a diagonal version with fourth order implicit 
smoothing terms as explained in chapter III. The flow solver converges fast to about four 
orders of magnitude reduction in RMSDQ. The computing time per iteration per grid 
point is 91 fxstc for the block tridiagonal version and 32 fx sec for the diagonal version of 
the code. 

In figure 5.16, particle traces for a single post at Re= 1000 are shown. A saddle 
point separation and a horseshoe vortex can be seen from the traces near the flat plate. 
The secondary flow in front of the cylinder wraps around toward the wake region and 
forms a counter-rotating pair of vortex filaments. These spiraling twin vortices demon- 
strate a striking difference between this type of juncture flow and a 2-D cylinder. The 
vortex filaments are washed upward. They attenuate as they interact and move down- 
stream. In reality, vortex shedding and possible unsteady motion take place at this stage. 
These tornado shaped vortices are very difficult to observe experimentally, and validation 
of this phenomenon was very much needed. Recently, Schewe (Schewe, G., Private com- 
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munication, 1985, DFVLR, West Germany) produced an oil flow visualization around a 
single post which, shows a clear evidence of the twin vortex behind the cylinder as shown in 
figure 5.17. This experimental observation is qualitatively similar to the computed results, 
shown by the particle traces in Figure 5.16. This juncture flow structure will lead to a 
strong variation in skin friction and pressure along the cylinder and hence significantly 
affects the overall loading on the post. 

5.0.2 Wing-Flat Plate 

Another si mil ax juncture flow, namely the wing on a flat plate case, was studied by 
Burke (1989). Numerical results were obtained by applying INS3D and compared with 
the experimental data of Dickinson (1986). The wing shape is a hybrid consisting of a 
1.5:1 elliptic nose and a NACA 0020 tail joined at the location of maximum thickness. 
This is a generic configuration characterizing a wing-fuselage juncture of aircraft or a hull- 
appendage juncture of ships. The Reynolds number of the experiments, based on the chord 
length of the wing, is 5.0 x 10 5 . The coordinate system is shown in figure 5.18. 

Turbulence modeling for this type of complex flow is very challenging. A high level 
turbulence model may be necessary for detailed study on the juncture flow itself. Since 
this study was done to extend the flow solver to more realistic applications, computing 
efficiency is of major importance. Therefore, the computational simulation was done by 
devising a simple algebraic turbulence model derived from Patankar et al. (1979). The 
results compares remarkably well with experimented data capturing all important features 
of the flow. Some of the results are reproduced here. In figure 5.19, surface pressure on 
a flat plate near the wing is compared with measurement while, in figures 5.20 and 5.21, 
velocity contours at two different vertical planes axe shown. Overall, the computed results 
compare very favorably with the experimental data. 


So fax several representative validation cases have been presented in this chapter. 
Other validation cases can be found in the references cited throughout this note. In the 
next chapter, two important recent applications of these codes will be presented. 
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VI. APPLICATIONS 


As discussed in Chapter I, there are a wide variety of 3-D problems where the 
present flow solvers can be utilized to analyze current configurations and to verify new 
design concepts. In this chapter, computed results of two major 3-D applications will be 
discussed. 


6.1 SSME Power Head Flow Simulation 

As the role of the Space Shuttle becomes increasingly important in scientific and 
commercial applications, it is desirable to increase the payload capability. Recently, an 
upgrade of the SSME power head was initiated to substantially increase the operating 
margin and the engine durability. To achieve this goal without increasing the weight and 
size of the existing engine, it became essential to understand the dynamics of the hot-gas 
flow in the engine power head. Because of the complexity of the geometry, an experimen- 
tal approach is extremely difficult as well as time consuming. Computational simulation, 
therefore, offers an economical alternative to complement the experimental work in an- 
alyzing the current configuration, and to suggest new, improved design possibilities. In 
the past few years, major milestones have been established from the computational effort. 
Here, highlights of our initial task are presented. 

6.1.1 Background of SSME flow analysis 

In the SSME staged combustion cycle, the fuel is partially burned at very high 
pressure and relatively high temperature in the preburners. The resulting hot gas is used 
to run the turbine and is then routed to the main injector where, along with additional 
oxidizer, it is injected into the main combustion chamber. The Reynolds number of the 
primary flow in the manifold is of the order of 10® per inch. Because of the high gas 
temperature, the Mach number is less than 0.12. The flow is turbulent and is practically 
incompressible. 

Figure 6.1 illustrates the current arrangement for the SSME power head components. 
Hot gas discharged from the gas turbine enters the annular turnaround duct (TAD) and 
experiences a 180° turn before it diffuses into the fuel bowl. This assembly is called the 
hot gas manifold (HGM). The gas flows into the main injector through three transfer ducts 
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on the left side of the power head (fuel preburner side) and enters into the region of the 
main injector posts. On the right side of the power head (oxidizer preburner side), there 
are two transfer ducts connected to the right side of the main injector assembly. Around 
the bottom portion of each liquid oxygen (LOX) post in the main injector assembly, there 
are a number of small holes through which the hot gas flows into the main combustion 
chamber. There it mixes with the oxidizer, which comes through a circular passage along 
the centerline of each LOX post. As a part of the engine development effort, a CFD study 
has been conducted to simulate the dynamics of the hot-gas flow in the power head. At 
the time this effort was started only the first version of INS3D was available. Therefore, 
the computed results presented in this section have been obtained using INS3D code. 

6.1.2 The computer model and the grid 

A computational model of the power head is chosen to analyze critical areas where 
the dynamics of the hot-gas flow are expected to have a significant effect on the overall 
performance of the SSME. As shown in figure 6.2, the model starts from the g as turbine 
exit on the fuel preburner side, and extends to the main injector assembly. The main 
injector consists primarily of a bundle of LOX posts, which is physically modeled by a 
porous media. 

Figures 6.2(a) and 6.2(b) demonstrate the 3-D grid for the SSME HGM. They show 
a horizontal and a vertical cross section, respectively, of the HGM. Figure 6.2(c) illustrates 
the details of H-grids for the cross-section of the three transfer ducts. This H-grid is 
generated for a unit circle. Near the boundary the grid lines are concentric circles except 
in the vicinity of the four singular points. Using the nearly orthogonal grid in this unit 
circle, one can obtain H-grids for tubes or ducts of any given shape and dimension by a 
simple linear transformation. Figure 6.2d shows an unwrapped surface of the annular fuel 
bowl with openings. The elliptical shape shown in figure 6.2c represents a cross section 
of the transfer duct of an advanced two-duct design, which will be explained later in this 
section. The current H-type grid topology for the circular cross-section was chosen to 
allow smooth transition from the axisymmetric TAD to the transfer ducts. Also, this 
arrangement is convenient in clustering the grid near the duct wall. 

The grid for the entire HGM system is generated by using algebraic functions, and 
is written with a high degree of flexibility for changing geometric configurations. By 
specifying the shape, the dimension, and the desired number of transfer ducts, a grid for a 
variety of new HGM configurations can be obtained in a short time. The ducts described 
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in this note axe connected directly to the fuel bowl without any fairings, while in the 
current engine the three transfer ducts are connected smoothly to the annular fuel bowl 
with fairings. This configuration with an abrupt change in geometry is more demanding 
computationally than smooth configurations. 

0.1.3 Multiple-Zone Computation 

A large number of mesh points are required to solve the 3-D turbulent flow in the 
SSME. To facilitate numerical simulations, the domain of interest is divided into several 
zones. This requires a special treatment at the interface for a smooth continuation of 
the solution between zones. Figure 6.3 illustrates a five-zone arrangement for the HGM 
flow field. Zone 1 includes the TAD and fuel bowl. Zones 2, 3, and 4 comprise the three 
transfer ducts for the current engine configuration. The annular region denoted as the 
racetrack of the main injector is represented by Zone 5. Also shown in the figure are some 
overlapping grids in the various zonal interfaces. The grid is chosen to be continuous and 
smooth across the zonal boundaries. In the first attempt at simulating this model, the 
racetrack (Zone 5) was not included in the computation. In a later computation involving 
a new design configuration, the racetrack region was included as well as the main injector 
assembly which was modeled by a porous medium. More detailed study of the LOX post in 
the main injector assembly was also performed using the results of the HGM computation 
as the boundary condition but will not be given here. Since the vertical plane through the 
center of the fuel bowl and the main injector is taken to be a plane of symmetry, results 
for the HGM presented in this note were obtained using only half of the HGM. When the 
swirling flow from the gas turbine exit is included, the symmetry assumption should not 
be used. 

In the pseudocompressible formulation, waves are propagating in both up- and down- 
stream directions while the solution approaches a steady state. In the present problem, the 
interfaces between zones are locations where the geometry changes abruptly. Therefore, in 
the neighborhood of those interfaces, the flow is expected to experience a rapid change. To 
maintain a smooth continuation of the solutions across these zones, and hence to achieve 
a stable and fast-converging computation, a means for providing adequate communica- 
tion for the traveling waves must be established. Overlapping regions and a proper zonal 
interpolation scheme are thus required for this purpose. 

A forward or backward differencing, if applied to the interfaces of multiple zones, 
would distort the geometric representation. To maintain a smooth transition of the flow 
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field across a zonal boundary, the Jacobian and the metrics at the interfaces are computed 
using grid points in neighboring zones. Then the pressure and the velocities, Q, are updated 
explicitly at each iteration. Let values at (n + 1/2) denote conditions to be used to advance 
the computation to n + 1. The values of Q n + l l 2 for Zone 1 at the exit plane are obtained 
from the values of the corresponding plane of Zone 2 at n, i.e., 


Q 


n+l/2 

B.C. 


Zone 1 


[Q interior] Zone 2 


Values of Q n+1 / 2 for Zone 2 at the zonal interfaces are taken from the latest computed 
result of Zone 1 as 

When more than two points are overlapped, the latest values in the interior of this over- 
lapping region must be properly transmitted to the next zone. There are a number of ways 
to treat this problem. The simplest one is take an average of the two values computed in 
Zones 1 and 2, as below: 

'«" +,/2 ] = \ ([«■]«. + W" + 1 W.) 


A scheme using updated zonal boundary values, but using original interior values, 
has also been tested. Either way, converged steady-state solutions have been obtained. 
However, the scheme with interior updates converges at a much faster rate. 

6.1.4 Grid effect 

The present flow solver (INS3D code) has been validated by computing fundamental 
fluid dynamics problems as illustrated in the last chapter. There are, however, many other 
aspects to be clarified in the real world applications, such as the grid induced error due to 
skewness and 3-D stretching, and grid clustering and its relation to numerical dissipation. 
In retd geometry, the flow is likely to go through strongly curved sections. This will 
introduce yet smother subject, that of strong curvature effect on turbulence structure. 

In the present application, for continuity and smoothness across zonal boundaries, 
an H-type grid is chosen for the circular transfer duct as shown in figure 6.2c. This is then 
connected to the side- wall grid of the HGM sis illustrated in figure 6.2d. In generating this 
grid for the SSME, errors are introduced mainly due to grid singularities, skewness, and 
stretching. 
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The present AF algorithm integrates the difference equations along the transformed 
coordinates. At the junction of the two H-grid directions, flow particles in the two co- 
ordinate directions could communicate only indirectly via the interior mesh points. This 
produces some corner-effect error. Even though this error is not as severe as the one caused 
by external flows, an ad hoc method of eliminating this corner effect is devised based on 
a finite-element concept. Let j and k denote the indices for the grid points along the 
increasing £ and 77 directions, respectively. Then let j = k = 1 be the corner point, which 
is singular. First, the pressure at this point, pn, is determined by an extrapolation along 
the diagonal direction j = k. Second, p\$ and P31 are obtained in the usual manner. Then 
P12 and P21 are established by an interpolation along the circular surface. 

The full viscous term given by equation ( 3 . 18 a) can be simplified to equation ( 3 . 18 b) 
when the grid is orthogonal. Even though full viscous terms can be used, it is convenient 
and economical to keep only the orthogonal part. It is, therefore, of practical interest to 
estimate the overall error caused by the orthogonal formulation when a nonorthogonal grid 
system is used. 

As a quick measure of an overall error caused by grid skewness, a 2 -D channel flow 
is computed using two grids, namely, 1) a stretched Cartesian grid (orthogonal) and 2) 
a nonorthogonal grid where the skewness is controlled on the upper half of the channel 
as sketched in figure 6 . 4 . In the computation only orthogonal terms are kept. Converged 
solutions on the lower half, where the two grids are identical^ are then compared. Total 
error depends additionally on the Reynolds number and the third- directional skewness. 
However, this quick experiment indicates that the orthogonal assumption can be used 
without significantly impacting on the overall solutions if the grid is basically orthogonal. 
Full investigation should include more severe cases involving separation and recirculating 
regions. 

0.1.5 Turbulence models 

The quality of the solution under turbulent conditions will be validated next. The 
flow through the SSME power head offers a variety of rich internal flow phenomena. Among 
others, the curvature effect of turbulence was most uncertain at the time the current 
simulation started. Therfore, the validation effort which focused on this will be summarized 
here. 


The current configuration of the SSME hot gas manifold was designed including 
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an axisymmetric annular 180° TAD. In the sharp U-turn region, the ratio of the internal 
shear layer thickness to the radius of curvature is of order one. On the convex side, i.e., 
the inner wall, the turbulence is expected to be greatly reduced, while on the concave side, 
it is expected to be subst anti ally enhanced. Furthermore, strong adverse and favorable 
pressure gradients coexist and interact with each other. Because of the sharp turn, the 
streamwise variations and the normal gradient of the streamwise velocity are of the same 
order of magnitude. Furthermore, the flow rapidly changes direction. Understanding the 
structure of turbulence in this flow is crucial for a successful computation of the TAD. 

Since the original work by Prandtl in 1929, the effect of mild steamline curvature 
on the structure of a turbulent boundary layer has been studied by many investigators 
(for a comprehensive review, see, for example, Bradshaw, 1973). It is well known that a 
mild curvature produces a profound effect on the turbulent flow structure. For example, 
if the ratio of the boundary layer thickness to the radius of the wall curvature is merely 1 
percent, there will be approximately 10 percent or more change in the turbulent quantities. 
In most external flows such as the flow over an airfoil, the streamline curvature is designed 
to be very mild except near the nose region where flow remains laminar. For problems 
of this kind, the literature on the curvature effect is quite extensive (see, for example, 
Wattendorf, 1935; Eskinazi and Yeh, 1954; Gillis and Johnston, 1983; Moser and Moin, 
1984). However, study on a strong curvature effect is very limited. 

To meet this demand, several experimental investigations were started to study 
the fundamental structure of turbulence associated with sharp U-turn internal flows. In 
one experiment, which was performed at Rocketdyne by Sharma et al. (1987), turbulent 
flows in an axisymmetric U-duct were studied. Figure 6.5 shows the TAD model for this 
experiment. The duct width is 2 inches, and the radius of the centerplane of the annular 
duct is 10 inches upstream of the turn and 14 inches downstream of the 180° bend. The 
radius of curvature is 1 inch for the convex inner wall, and 3 inches for the concave outer 
wall. From the beginning through the completion of the turn, the annular cross-sectional 
area undergoes a 40 percent increase. The boundary layer on the inner wall side is therefore 
expected to separate in the neighborhood of the 180° turn. The experiment is conducted at 
atmospheric pressure. The Reynolds number is approximately 10 s , and the Mach number 
is about 0.1. Hot wire and hot film probes are used for data acquisition. 

For such complicated turbulent flow as encountered in the SSME, high level turbu- 
lence models such as a two-equation k-e model or Reynolds-stress model may be needed. 
However, for many engineering applications, a simple yet adequately accurate model will 
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be of considerable value for design purposes. Therefore, developing and validating a simple 
algebraic turbulence model which can be used for simulating flows with strong streamline 
curvatures is currently of interest. This case is used as a part of the general flow solver 
and turbulence model validation for application to the SSME. 

Several levels of turbulence models have been implemented into INS3D. These in- 
clude the Baldwin and Lomax (1978) algebraic model in which length scale is determined 
by the location of the maximum moment of vorticity. This model has been widely applied 
in the external flow problems. However, the maximum moment of vorticity is not as well 
defined for fully turbulent internal flows as for external flows. In particular, the moment 
of vorticity is almost constant for a fully developed pipe or channel flow except in the 
sublayer region. 

For fully turbulent internal flow in a duct, the two boundary layers associated with 
the opposite walls merge together. For these flows, there exists a pair of opposite- sense 
vortices. Therefore, there exists a position where these two opposite vortices cancel each 
other. This position is clearly the one which divides the two boundary layers. Chang et al. 
(1985b) proposed a simple extended Prandtl-Karman mixing length theory which, together 
with the strength of vorticity, forms an eddy viscosity model. The vorticity length scale is 
defined to be the distance between the location of maximum vorticity and the point where 
vorticity first vanishes. This length is an active measure of the turbulent eddy size. The 
location of the maximum vorticity for an attached boundary layer is at the wall. When 
the layer begins to separate in the presence of a strong adverse pressure gradient, the 
point where the magnitude of vorticity is maximum starts to move away from the wall. 
The point of vanishing vorticity can be replaced by the location of minimum vorticity in 
merged layers for the purpose of a numerical study. The accuracy of this model is evaluated 
here using the present test case. 

In the present experimental setup, it is very difficult to produce a fully turbulent flow 
in the entire duct. Even with a boundary layer trip mechanism, the artificially thickened 
layers upstream of the turn cover only about half the duct width. In the center core region, 
the flow remains practically inviscid. This inviscid core flow acts like a free-stream in an 
external flow. As free-stream conditions are approached near the edge of the boundary 
layer, an intermittent zone of turbulent and nonturbulent flow is encountered. This and the 
wake component in the outer layer change the structure of the shear layer drastically from 
that of a fully turbulent flow. For these flows, the mixing length is modified accordingly. 
For more details on the model, see Chang and Kwak (1988). A review on various levels of 
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turbulence models can be found in articles by Reynolds (1976) or Ferziger (1987). 


In the present problem, only three azimuthal planes axe required because of the 
axisymraetry. The center plane is at 8 =0 and other two planes are set on either side 
of the center plane offset by a small angle. Across the duct width, 81 grid points are 
distributed with an expansion factor of about 1.2 . The smallest grid size next to the wall 
is 0.02% of the duct width which corresponds to y + of less than 5. In order to resolve the 
detail of the separated flow profile, the largest grid size in the core region is set to be 0.025 
times the width. 

The 180° U-turn region is divided by 42 equal intervals. An orthogonal cylindrical 
grid system is used here to eliminate any uncertainty in accuracy due to non-orthogonality. 
The grids upstream and downstream of the bend axe then fixed by smooth expansions from 
the bend region. Close to both walls these grids are orthogonal. The flow path is therefore, 
represented by a 3x81x95 grid in the circumferential, radial and flow directions respectively. 
This grid in the U-tum region is shown in figure 6.6. 

Figure 6.7 shows the comparison of the static pressure distributions between the 
measured data and the computed results. The measured data were obtained in an earlier 
experiment in which no turbulent tripping mechanism was used. The resulting boundary 
layer on the outer wall at 8 inches upstream of the bend, i.e., x/h = -4, is very thin. 
By using these measured velocities as an inlet condition at that location, the computed 
result shown with a solid line agrees very well with the experiment. In the neighborhood 
of the 180° turn, pressure on the inner wall shows a dip followed by a long recovery from 
the trough. In external flows, boundary layer separations usually interrupt the pressure 
recovering processes. The resulting pressure downstream of the separation will be almost 
constant. In the present confined duct, the flow which is pushed outward due to inner layer 
separations is squeezed back by the flow on the opposite side. As a result, the effective 
flow area is reduced creating a region of local acceleration which results in a local pressure 
dip. 


In a more recent experiment, the layer is tripped. It is artificially thickened to 
about 0.5 inches at the position 10 inches upstream. Because of the thickened layer at the 
outer wall, a larger portion of kinetic energy is carried by the inner-wall region inducing a 
larger acceleration in the bend. As a result, the static pressure on the inner convex wall is 
expected to be a little lower than in the previous run. 
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Comparisons of the measured and the computed streamwise velocity profiles are 
shown in figure 6.8. Agreement is in general very good except near the outer concave wall 
in the bend region. The numerical results show a slight separation from the beginning of 
the turn through about 103°, while the experimental data shows attached layers all the 
way through the turn. Accuracy of the measured velocities in this region is still uncertain 
due to two factors. First, because of the sharp concave wall effect, the enhanced turbulent 
intensities obtained are as high as 45 percent of the measured mean velocities. Since the 
hot-wire probe measures only the rate of heat transfer, it does not detect the direction of the 
flow. In a region where heat transfer is significantly affected by the turbulent fluctuations, 
it is extremely difficult to obtain the true mean flow velocities as evidenced by experiments 
in the low speed region of a mixing layer. Second, the mass flow obtained by integration of 
the measured velocities over the cross-sectional surface of these few planes is substantially 
larger than the amount of mass coming from the inlet. The computed velocity profiles 
conserve the mass flow to an order of 10 -5 . Because of these uncertainties, no conclusion 
can be drawn for the discrepancies at this time. On the convex inner-wall side, the layer 
begins to separate at about 125° into the turn. The data obtained by the hot-wire probe 
registered a high value due to large turbulent fluctuations in the separated region, and the 
data do not reveal the direction of the flow. 

Overall, the accuracy of the solution using the present algebraic turbulence model is 
quite acceptable. Therefore, this model has been used quite extensively for simulating the 
SSME power head. In the remainder of this section, flow solutions in a variety of different 
HGM configurations at various Reynolds numbers are presented. Here, the Reynolds 
number is based on the mean velocity and the duct width at the entrance of the TAD. 

6.1.0 Current Three-circular-duct HGM Analysis 

First, steady-state solutions are obtained for the current three-circular-duct HGM 
at Re=1000. As illustrated in figure 6.9, for the initial computation, the racetrack of the 
main injector is not connected to the ducts. The three transfer ducts are assumed to 
discharge the flow separately, which results in no communication of the pressure between 
the center and the outer ducts at their exit planes. Therefore, the downstream condition 
is not a good approximation to the real case. For this reason, small residual waves have 
remained in the computed results. However, the root-mean-square value of the change 
in the flow variables, RMSDQ, has dropped below 10 -5 , and an essentially steady-state 
solution has been obtained. 
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In figures 6.10(a) and 6.10(b), velocity vectors are shown in the horizontal and 
vertical cross-sections corresponding to figure 6.2. The flow in the center transfer duct is 
highly nonuniform, and a large separation region is formed just downstream of the entrance 
to the transfer ducts. By comparison of the vector length in figure 6.10(a), the flow in the 
center duct is much slower than it is in the outer ones. The results shown in these figures 
agree qualitatively with the airflow test data conducted at a Reynolds number of about 
10 6 . The predicted mass flow through the center duct is 9.8% of the total mass flow, which 
agrees with the test data. 

Figures 6.11(a) illustrates the 3-D velocity vectors at an unwrapped plane near 
the inner wall of the fuel bowl. A reverse-flow pattern is clearly visible near the inner 
wall. Three-dimensional swirl patterns are predicted in the vicinity of the entrance to 
the transfer ducts. Figure 6.11(b) is a photograph that indicates, by means of surface- 
streak (shear-pattern) visualization, the similar swirls at the corresponding locations in 
the airflow test. 

The existence of the swirls can be explained as follows. The flow coming from below 
has a large momentum due to the relatively small width of the annular duct. Among the 
streamlines of this flow in between the two ducts, there exists a dividing streamline. This 
streamline has a stagnation point at the top of the fuel bowl as shown in Figure 6.11a. 
On the left side of this dividing streamline, the flow is bent leftward to the center duct. 
Because of symmetry, a rightward flow is also approaching the center from the other side. 
When these opposite currents approach each other, another dividing streamline is formed 
with a stagnation point, again at the top wall. The stagnation pressure forces the streams 
to bend downward, and at the same time, the streams make a right-angle turn into the 
circular duct. Conservation of momentum thus requires the formation of swirls. 

The pattern of the swirl and its center depends on the relative strength of the ap- 
poaching currents. Near the center duct, double swirls of equal strength are formed because 
of symmetry. In the vicinity of the entrance to the outer duct, the current approaching 
leftward from the rear part of the bowl is larger than the one approaching rightward. A 
stronger swirl is thus formed, locat.ed sideways toward the weaker stream. 

Figures 6.12 and 6.13 are the perpendicular cross-sectional views showing three 
different sections of the transfer ducts; namely, near the entrance, at the midsection, and 
near the exit plane. Near the entrance (figures 6.12a and 13a), the velocity vectors in 
the center duct have symmetric double swirls, while the outer duct has a strong swirl 
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accompanied by a much weaker one. The swirling velocities are largely reduced at the 
midsection and are physically dissipated before entering the main injector regions. 

0.1.7 Development of New Design Configuration 

From this computational flow analysis and also from experiments, the center duct 
of the current three-duct HGM is found to transfer a limited amount of mass flow (about 
10% of the total flow). Also the transverse pressure gradient remains large with a large 
bubble of separation after the 180° turn. To improve the quality of the flow, a large-area, 
two-duct design concept has been developed. In addition, the ducts are chosen to have an 
elliptical shape in order to distribute the mass flow evenly to the main injector region. An 
outer-surface grid for a two-duct model is illustrated in figure 6.14. 

First, to reduce the large separation bubble after the 180° turn, a parametric study 
is performed to find the best possible configuration. In figures 6.15(a) and 6.15(b), com- 
parison of the current three-duct configuration and the new two-duct design is illustrated. 
As shown in figure 6.15(a), a large separation bubble existing in the present design is prac- 
tically removed in the new configuration, shown in figure 6.15(b). This is confirmed by 
experiment as shown in figures 6.15(c) and 6.15(d) for the current and new designs, respec- 
tively, where velocity measurements at five different locations across the channel between 
the inner and the outer wall are shown. To find the most favorable flow conditions, over 20 
different two duct configurations were studied computationally, thus providing potentially 
the optimum geometry to the designers. Ideas from researchers were incorporated into the 
current human optimization process. Fully numerical optimization may be considered in 
the future in conjunction with development of a faster code. 

Figure 6.16(a) illustrates an example of the triple-swirl pattern in the elliptical duct 
near the entrance at Re=10 3 . Here, because of the absence of the center duct, the stream 
approaching from the left side is much stronger than in the previous case. Another point 
of interest is that the upcoming stream entering the elliptical duct directly from below 
is also more massive. The three currents tire of almost the same strength, resulting in a 
triple-swirl flow. The swirling is greatly dissipated along the duct. Figure 6.16b shows the 
remaining small swirling vectors at the duct exit. A steady-state, turbulent-flow solution 
for an HGM with two elliptical transfer ducts at Re = 10 5 has been obtained by Chang 
et al. (1985). In figure 6.17, the turbulent swirling flow in the elliptical transfer duct is 
illustrated. In comparing these results with the l amin ar solution in figure 6.16, it is seen 
that only a double-swirl pattern exists. 
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Once the new configuration is developed, it is necessary to analyze the new engine 
relative to the old one. To simulate the new two-duct engine accurately, a more complete 
numerical model was developed as shown in figure 6.18 where both the fuel side (left side) 
and the oxidizer side (right side) of the power head are included in the computer model. 

One important objective of the redesign was to reduce the transverse pressure gra- 
dient at the exit of the gas turbine underneath the preburner. Therefore, the turbulent 
solution for the new two-duct HGM is compared to that of experiments performed using 
both the current three-duct HGM and the new two-duct configuration. As shown in Figure 
6.19, the pressure gradient around the fuel bowl of the HGM is greatly reduced from that 
of the original configuration. From the hardware perspective this results in substantially 
reduced load on the bearings which hold the gas turbine and the turbopump sketched in 
figure 6.1. 


The most significant objective of the present study is to pinpoint the locations 
where flow experiences the greatest energy losses. An important measure of the energy 
losses is the mass-weighted average total pressure along the flow. Figure 6.20 illustrates 
the decreasing coefficient of the mass- weighted total pressure along the centerline of the 
TAD, the fuel bowl, and the transfer duct. The total pressure coefficient Cpo is defined as 

„ _ Po - Poi 
Cpo ~ ~PoT~ 


where 



dm 


The discontinuities shown in the figure correspond to the entrance of the duct where energy 
fluxes are computed over different planes. In the figure, three different HGM configurations 
are compared. The initial two-duct design shows 28% less total pressure drop compared to 
the current three-duct version. After fine-tuning the two-duct configuration computation- 
ally, the pressure drop decreased even further to 36% less than the original configuration. 
This final configuration is then tested using cold air flow, which shows 40% reduction in 
pressure loss. 


6.2 Artificial Heart Flow Simulation 

As shown thus far, the major advances in CFD have been made in aerospace en- 
gineering. With the advent of supercomputers as well as the developme nt of fast al- 
gorithms, computational flow simulations have become a practical means for aerospace 
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designs. Therefore, it will be of considerable benefit to medical researchers to extend this 
CFD technology to biofluid analysis. Some of the important studies in biofluid mechanics 
deal with aspects of blood flow, such as heart and blood vessel problems (see, for example, 
Peskin, 1982), which are relevant to cardiovascular diseases and their treatment. Under- 
standing the flow phenomena through numerical simulations will be of significant value 
toward finding treatments for these problems. Therefore, the potential payoff to human 
health will be tremendous. In this chapter, the potential of CFD in biofluid mechanics 
research is demonstrated by simulating the flow in an artificial heart. 

6.2.1 Background 

Recently, the demand has grown for mechanical hearts or ventricular-assist devices 
(VAD) for use as a temporary life support system. These devices are becoming a power- 
ful tool for assisting patients to recover from heart attacks or as a temporary bridge to 
transplants. Presently, There are several problems with these devices mainly involving the 
flow of blood. A research and development effort is in progress to improve the flow quality 
as well as to develop better material and control systems. Experimental investigations on 
these devices are very limited and many aspects of the flow are yet to be studied. There- 
fore, it will be very valuable to medical researchers to simulate the flow in these devices 
by applying state-of-the-art CFD technology. 

Blood flow through these mechanical devices is very complicated in many respects. 
The fluid may exhibit significant non-Newtonian characteristics locally and the geometry 
is usually very complicated. In an artificial organ, as red blood cells go through high- 
shear turbulence regions, they may be damaged; the downstream region of an artificial 
heart valve is an example. The flow is unsteady, possibly periodic, and very viscous and 
incompressible. This problem is very much interdisciplinary and an attempt for a complete 
simulation would be a very formidable task. However, an analysis based on a simplified 
model may provide much-needed physical insight into mechanics of the blood flow through 
these devices. 

* 

The formulation of the flow solvers described earlier in this note is based on a New- 
tonian fluid assumption. However, since the governing equations are solved in a generalized 
coordinate system, viscosity that varies in space and time as well as moving geometry is al- 
lowed. These flow solvers can be applied to analyze mechanical hearts and ventricular-assist 
devices. A full simulation of viscoelastic flow is very difficult because of the nonlinearities 
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of the fluid. However, as a first step toward full simulations, non-Newtonian effects of the 
blood flow can be simplified by a constitutive model for the viscous stresses. The primary 
purpose of the current application is to transfer technology developed for aerospace appli- 
cations to artificial heart reseaxch. This requires the state-of-the-art in CFD for treating 
unsteady interned flow with moving boundaries. Researchers in the aerospace field will 
benefit from the advances made in unsteady flow simulation while the artificial heart de- 
velopers benefit by gaining a better understanding of the fluid flow within their devices 
which hopefully will lead to an improved design. 

There are many variations of VAD and artificial hearts. Among those, two different 
approaches are illustrated in figures 6.21 and 6.22. In the device depicted in Figure 6.21, 
the blood is pumped by pressurized air supplied externally. Simulation of this device 
requires a moving boundary procedure where the boundary itself should be determined by 
the motion of two fluids. Another approach is developed by Pennsylvania State University, 
where the pumping is done by an electric motor. This heart is sketched in figure 6.22. Our 
demonstration calculation is performed on this Penn State artificial heart. Experiments 
on various configurations are in progress at Penn State University. In the computational 
simulation, the time-accurate INS3D-UP code is used. The geometry of a computer model 
and some preliminary solutions are presented in this section. 

6.2.2 Geometry and Grid Generation 

The actual model of the Penn State artificial heart poses some very difficult problems 
from a computational standpoint. A computer model of the heart, i.e., surface grid and 
a shaded-surface representation, is shown in figure 6.23. The heart is composed of a 
cylindrical chamber with two openings on the side for valves. The pumping action is 
provided by a piston surface which moves up and down inside the chamber. The actual 
heart has a cylindrical tube extending out of each of the valve openings. These tubes 
contain tilting flat disks which act as the valves. The current computational model will 
neglect the valves altogether and will use the right and left openings shown in figure 6.23 for 
the inflow and outflow boundaries, respectively. In the computations, as the piston reaches 
its topmost position, the outflow valve closes and the inflow valve will open instantaneously. 
Similarly, as the piston reaches its bottommost position, the outflow valve will open and 
the inflow valve will close. 

In the actual heart device the piston moves through the entire chamber volume, 
which includes most of the valve opening. This will cause some severe problems for the 
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computational grid as the piston moves across the valve boundaries. Since the flow solver 
is designed to use body-fitted coordinates, it is necessary to place the grid lines around the 
valve to coincide with the valve opening boundaries. Yet since the piston moves past this 
opening, the grid will have to accommodate both of these surfaces. There are several ways 
in which a computational grid can handle this motion. One method is with the use o.f a 
Chimera scheme (see Beneck et al ., 1985) in which two grids are used, one which moves 
with the piston, and one which is attached to the rest of the body. Information is passed 
between the two grids by interpolating variables at the grid interfaces. This method 
can be quite expensive and very complicated to implement. Two other methods which 
are somewhat simpler to implement require one grid inside the computational domain. 
One of these methods use a stationary grid through which the piston surface travels. 
Boundary conditions applied at the piston surface would ensure that mass and momentum 
are conserved for the partial grid cells as the piston moves through them. However, with the 
curved surfaces at the valve boundaries, this method would require partial grid cells with 
an arbitrary number of sides, for which generalizations are difficult. Another single grid 
method uses a grid which varies in time as the piston moves through the boundary. The 
generalized coordinate transformation enables provision for grid motion in the equations 
by use of the time varying metrics. This method requires that a separate grid be generated 
for each discrete position of the piston during the calculations. This may be the simplest 
method that will accommodate the entire piston motion, yet it does leave a difficult grid 
generation problem. 

In an effort to eliminate this problem, a simpler piston motion was decided upon 
for the calculations of the present study. The piston motion was restricted so that the 
piston did not rise above the bottom of the valve openings. The single grid approach that 
varies in time was used, so that a constant number of grid points are used, and the grid 
in between the piston and the bottom of the valve openings compresses and expands as 
the piston moves up and down. To facilitate this action, the piston was allowed to travel 
further down, so that the overall volume of the chamber in the computations was larger 
than in the actual device. 

To make the most efficient use of grid points, an H-H grid topology was used to 
fit the grid to the computational domain. The grid dimensions were chosen to be 39 x 
39 x 51. In order to generate a grid at each time step for this geometry, the surface grid 
was first generated, and from that, an algebraic grid generator and elliptic smoother was 
used to generate the interior points using the distribution given on the surface grid. To 
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generate the surface grid, the side boundary was divided into seven different zones. The 
grid points were distributed along each of the zonal boundaries, and then a biharmonic 
solver was used to generate the grid interior to each of the zones. The biharmonic solver 
was also used to generate an H-grid for the the top and bottom surfaces of the heart device. 
The inner points were then computed with the use of an algebraic solver coupled with an 
elliptic smoother. This approach made it relatively simple to repeat the process at each 
time step for any given position of the piston surface. 

6.2.3 Computed Results 

The calculations were carried out on a Cray 2 supercomputer with a core memory of 
256 million words. This large amount of memory has proven to be very useful because the 
implementation of the implicit scheme in INS3D-UP requires about 180 times the number 
of grid points in words of memory. For the current modest grid of 77,571 points, over 14 
million words are required. The computations started with fluid at rest, the piston at the 
bottom position, and the outflow valve open. The computations were carried out for a 
Reynolds number of 500. The pseudotime-step At was set to 10 12 and the physical time 
step At was set to 0.025. The velocity of the piston during the motion in between its 
maximum and minimum positions was set to be constant, which is very nearly the same 
as that in the actual heart device. The period of the entire piston cycle was set to a time 
of 5 nondimensional units, so that 200 physical time steps were required for each full cycle 
of the piston. 

During each time step, the subiterations were carried out until the maximum residual 
dropped below a specified small value or until a maximum of 20 subiterations were used. 
During most of the piston cycle only 12-15 subiterations were required, but when the 
piston was changing directions, it did not completely converge in 20 iterations. This did 
not cause any stability problems, yet it remains to be seen what effect this has on the 
accuracy of the solution. The computing time required for each period of piston motion 
was approximately four hours. The computations were run for four periods during which 
time particle paths were computed after being released near the inflow valve. 

Some very interesting flow physics were observed during this period of motion. Var- 
ious post-processing techniques have been utilized in analyzing the results, such as the 
velocity vectors colored by pressure level, and vorticity magnitude contours which can be 
related to the wall shear stress. Figure 6.24a shows particle traces as the piston nears its 
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bottom position. Two distinct vortices axe seen to have formed from the flow separating 
as it enters through the inflow valve. In figure 6.24(b), an experimental photograph (J.M. 
Tarbell, Pennsylvania State University, 1988) shows bubbles entering the inflow valve as 
the piston nears its bottom position. A similar two-vortex system is seen to form here. 
Figures 6.25a and 6.25b show velocity vectors during the inflow phase in planes passing 
through the center of the inflow valve. The first of these show a top view of vectors in a 
plane to the piston, while figure 6.25(b) shows a side view of vectors in a plane perpen- 
dicular to the piston. These figures portray the complexity of the vortical structure of 
this flow. Figure 6.25(a) again shows the presence of two vortices formed as the incoming 
flow forms a jet. Figure 6.25(b) shows also how the flow separates underneath the valve 
opening, although the flow there is strongly three-dimensional. Also in the figure is seen 
the presence of additional vortices against the back wall opposite the valve opening. 

Some measurements of turbulent flow through the artificial heart geometry were 
taken by Tarbell et al. (1986). They were able to measure some of the mean flow charac- 
teristics, as well as turbulent shear stresses. However, these measurements were taken with 
the disk valves in operation. These valves greatly affect the incoming fluid by directing 
it along the wadis so as to set up one large vortex inside the chamber which washes the 
walls free of any stagnating fluid. This effect is by design so that the blood would not tend 
to clot due to a large residence time in a stagnation region. Since these valves were not 
included in the current calculation, no qualitative comparison can be drawn between the 
experiment and the current calculation. 

The present solution shows the capability of the computational procedure for sim- 
ulating complicated internal flows with moving boundaries. For this initial calculation, 
the motion of the piston of the actual device was simplified. Also simplified was the valve 
opening and closing procedure. One difficulty in performing the present study is that there 
is currently very little validation available for these computations. It is realized, however, 
that this work represents a first step toward developing a CFD tool for this type of flow. 
With the use of more advanced grid techniques, and multiple zones to handle the valves, 
results of the entire artificial heart project will be reported in the future. 
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VII. CONCLUDING REMARKS 


In the present note, numerical simulation methods for viscous incompressible flows 
are discussed from an application point of view. Our main interest has been in the 3-D 
real-world geometry. Naturally, the computational requirements for these problems are 
different from the fluid dynamics studies in the fundamental category. Therefore, the 
main emphasis has been placed on a primitive variable formulation which has been the 
most commonly used approach for 3-D problems to date. Detailed discussion has been 
given on flow solvers developed at NASA Ames Research Center. The validation cases and 
applications presented are the results from these solvers. 

There are wide spectrum of algorithms and procedures for incompressible flows not 
covered here. The number of numerical studies of various incompressible flow problems are 
also staggering. It is an almost impossible task to review all this material. The material 
presented here is that familiar to the author; thus many important results and references 
are naturally left out, some for brevity and some due to ignorance. Many comments are 
made based on the author’s experience without substantiation with firm numerical proof; 
they do not represent a consensus. 

There Me several important aspects left out in the discussion. In differencing the 
governing equations, conservative properties have been limited to momentum and mass 
conservation. Another property to be considered is the kinetic energy conservation. This 
is especially important in simulating turbulence involving small scale motions. In general, 
for the type of applications we have been interested in, the total number of mesh points 
and the total computing time limitations did not allow very fine mesh computations. 
For the class of problems requiring a fine mesh, such as the large eddy simulation of 
turbulence, kinetic energy conservation becomes crucial so a s not to numerically create 
small-scale turbulence, which will eventually lead to nonphysical accumulation of energy 
in high-frequency motion. 

In a more practical sense, even though computer speed and memory have been 
increased substantially in the recent past, the speed and the memory requirements of a 
flow solver are still dictated by the turn around time. In many engineering applications, 
it is very important to generate solutions in a timely fashion to have any impact on the 
design and analysis. We feel that numerical simulations can now provide complementary 
information to measured data, thus reducing the number of experimental trials required 
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for developing advanced flow devices. A more sophisticated approach, such as computer 
optimization of a design, can be attempted when the flow solver becomes at least one 
order of magnitude faster, or the computer speed improves that much perhaps by massive 
parallelism. 

There are several areas to be considered in developing a new generation of fast flow 
solvers. Several additions to the technology presented here will help to make the current 
methods as efficient as possible. Multigrid acceleration is one possibility which is physically 
consistent with the incompressible formulation. Use of a solution adaptive grid is also a 
practical approach to improve the accuracy and efficiency of the computation. This, as 
well as problems involving a moving domain, may require the generation of a new grid at 
each time step. For generating a grid at each time step, automated grid generation would 
be ideal. Overall, the solution procedure should be developed to best utilize computer 
characteristics such as vectorization, parallel processing, and access to memory. * 

Regarding the development of a universal code, it is hard to devise the best scheme 
for all flow speeds and for many different types of flow. Considering the uncertainties 
related to the turbulence model, it seems better to have a specialized code suitable for 
each different class of problems. For example, one can imagine having a code for an 
impeller type problem while having another for flow over aerodynamic shapes. Despite 
the limitations in algorithm speed and accuracy, and computer speed and memory, when 
the state-of-the-art flow solvers are utilized by creative researchers, the result will be 
tremendously beneficial for developing modern devices requiring viscous incompressible 
flow analysis. The current work in INS at NASA-Ames will continue to develop each of 
the codes presented in this note, which are considered to be the most promising methods 
for real-world flow simulations. As these codes are applied to these problems, we hope to 
quantify the advantages and shortcomings of these codes. 
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Figure 5.1 Channel flow : velocity profile at Re=1000, /3 —5, r— 0.1 
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Figure 5.2 Pressure contours for developing channel flow: (a) (3=5, (b) (3—0.1 



Figure 5.3 Convergence history for channel flow at Re— 1000: 
(a) RMSDQ, (b) RMSDIV 
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Figure 5.4 Geometry of backward-facing step flow problem. 
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Figure 5.5 Separation and reattachment 
lengths for the flow over a backward- 
facing step. 
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Figure 5.6 Convergence history for the flow 
over a backward-facing step : (a) residual, 
(b) primary reattachment length. 
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Figure 5.7 Time evolution of separation length for flow over 
a circular cylinder at Re=40. 
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Figure 5.8 Vortex shedding from a circular cylinder. 
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Figure 5.9 Karman vortex street behind a circular cylinder at Re=105: (a) computed (INS3D-UP), 
(b) experiment (S. Taneda, see Van Dyke, 1982) 


Figure 5.10 Geometry of a square duct 
with 90° bend. 



Figure 5.11 Streamwise velocity profile for flow 
through a 90° bend: (a) xy plane 
at z=0.25, (b) xy plane at z=0.5. 
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Figure 5.13 Predicted radial flow in 90° bend (McConnaughey et al., 1989, using 
INS3D) compared with the data of Taylor et al. (1982). 
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Figure 5.16 Particle traces for a single post at Re=1000. 



Figure 5.17 Oil flow visualization around a single post at Re=l. 85x10 s 
(by G. Schewe, DFVLR, 1985). 
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Figure 5.20 Velocity contours for wing-flat plate at x/c=0.18 and Re=5xl0 5 : 

(a) experiment (Dickinson, 1986), (b) computation (Burke, 1989; INS3D) 



Figure 5.21 Velocity contours for wing-flat plate at x/c=0.75 and Re=5xl0 5 : 

(a) experiment (Dickinson, 1986), (b) computation (Burke, 1989; INS3D) 
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Figure 6.1. SSME power head component arrangement 



Figure 6.2. Grid of the SSME hot-gas manifold. 

a) Horizontal view (cross section A-A); b) vertical view (cross 
section B-B); c) H-grid for circular and elliptic cross section 
of transfer duct; d) unwrapped surface of annular fuel bowl (cross 
section C-C). 
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Figure 6,3. Multiple-zone arrangement 
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Figure 6.4. Grid effect on channel flow solution at Re=1000. 

a) Definition of grid skewness; b) Relative error due to skewness. 
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(b) 


Figure 6.11 Velocity vectors on unwrapped surfaces: 
(a) unwrapped inner wall, (b) experiment. 
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Figure 6.12 Velocity vectors at cross 
sections of center duct in three-duct 
Hot Gas Manifold. 


Figure 6.13 Velocity vectors at cross 
sections of outer duct in three-duct 
Hot Gas Manifold. 
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Figure 6.14 Outer-surface grid for a two-duct HGM 
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Figure 6.15 Turn-Around-Duct redesign: 

a) Computed velocity vectors for the current design; 

b) Computed velocity vectors for the new design; 

c) Experimental velocity head measurements for the current design 
(6 defined in Figure 6.2a); 

d) Experimental velocity head measurements for the new design. 
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Figure 6.16 Velocity vectors at cross 
sections of transfer duct in two-duct 
HGM (laminar: Re=10 3 ). 


Figure 6.17 Velocity vectors at cross 
sections of transfer duct in two-duct 
HGM (turbulent::Re=10 5 ). 






Figure 6.18 Computer model of the new two-duct SSME power head : 

(a) vertical cross-section (B-B), (b) horizontal cross-section (A-A). 
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Figure 6.19 Pressure coefficient on fuel bowl outer surface after 180° turn. 
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Figure 6.20 Pressure losses in three-duct and two-duct HGM. 
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Blood out 



Figure 6.21 Schematic of an artificial heart operated by air pressure. 



Left diastole Left systole completed 

Figure 6.22 Schematic of the electric motor driven Penn State artificial heart. 
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(b) 


Figure 6.23 A computer model of the Penn State artificial heart: 

(a) surace grid, (b) main chamber showing valve openings. 
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a. Computational results 



lx Experimental results 


Figure 6.24 In comping particle traces as the piston nears the bottom position: 
(a) computation (Rogers et al., 1989), (b) experiment (J.M. Tarbell, 
Penn State Univ, private communication, 1988). 
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(b) 

Figure 6.25 Velocity vectors of incoming fluid: (a) top view in horizontal 

plane through center of inflow valve (cross-section A-A), (b) side 
view in vertical plane through center of inflow valve (cross-section B-B). 
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